Compiling the data

For this project we aimed to compile data on a large set of species in the Passerine bird order that would allow us to account for many evolutionary and environmental predictors of extinction. The following details the source and use of the data.

Speciation Rate

To obtain speciation rate estimates we used Bayesian Analysis of Macroevolutionary Mixtures (BAMM). Additionally, we compared our BAMM output files (event data) against a pre-existing extinction rate dataset (Harvey et al. 2017).

Rather than obtaining one estimate of speciation rate from one tree, we ran BAMM 100 times on 100 trees plus an MCC tree to obtain uncertainty estimates of the speciation rate generated from the phylogeny created by Jetz et al. (2012). We also estimated extinction rate from BAMM.

Sexual Selection

Our proxy for sexual selection was based on measures of sexual dichromatism. Notably, we use two existing datasets; a complete set of male and female plumage scores in Passerines from Dale et al. (2015) and reflectance data from 1000 birds from Armenta et al. (2008).

Briefly, Dale et al. (2015) obtained male and female plumage score, which can be divided to obtain a sexual dichromatism index (SDi). The plumage score was scored using the Handbook of the Birds of the World Del Hoyo et al. (n.d.). By scanning images of males and females in this book across multiple patches Dale et al. (2015) were able to obtain mean plumage scores from RGB values. Notably concern about the accuracy of these measurements were raised, thus this dataset was previously compared with a more thorough reflectance dataset of Australian birds, the two datasets correlated but was expectedly noisy. Here we compare the Dale et al. (2015) data against another reflectance dataset from Armenta et al. (2008) using colour discriminability as the parameter.

Alongside sexual dichromatism data there is also partial data for other metrics analysed in Dale et al. (2015) such as: Body size, tropical life history, Sexual selection, Cooperative breeding and Migration. We use a phylogenetic PCA (PC1) that estimates male-bias sexual selection from a species level of sexual size dimorphism, social polygyny and paternal care.

Environmental predictors

We obtained species range distribution maps from Birdlife International (through IUCN) (BirdLife International and Handbook of the Birds of the World 2017). These species range maps cover nearly all species of birds. From these spatial information files we were able to obtain estimates of the following (note that not all the data extracted was subsequently used in the analysis):

Range Size
Average and standard deviation in 19 bioclimatic variables (each range randomly sampled 1000 times)
Average and standard deviation in 19 bioclimatic variables from the last-glacial maximum (LGM) and the last-inter-glacial (LIG)
Net primary productivity (NPP) estimates and variability
Average and standard deviation of human population density in the species ranges.

The code for extracting the environmental data from the species is provided in another R-markdown document Github link, however no raw data is provided alongside this file due to the large file size of the shapefile and raster information. We have put the extracted environmental variables in a csv filed called complete.dataframe.csv for convenience. Briefly for each of the ~ 6,000 species, we:

1) Randomly sampled the range polygon 1000 times:

#Increse the iterations (defaukt is 4) so we can obtain complete samples of each range
bird.points <- lapply(bird.ranges, 
                       function(x) {spsample(x, n=1000, type="random", iter = 30)})

2) Extracted the bioclim or other (e.g. NPP) data from each of the points:

bird.values <- lapply(bird.points, function(x) {raster::extract(bioclim_data, x)})
bird.values <- as.list(data.frame((bird.values)))

3) Summarised the extracted data across the 1,000 points into a summary value of interest (means and se) for each species and exported that data.

#Obtain means
bird.frame <- lapply(bird.values, function(x) {as.data.frame(x)})
bird.summary <- lapply(bird.frame, function(x) {
(as.data.frame(apply(x,2,mean, na.rm =T )))})

#transpose
bird.means <- t(as.data.frame(bird.summary))
bird.means <- (split(bird.means, 1:19))

#Now add column of species name: Same order carried through
bird.means <- cbind.data.frame(shps.jetz, bird.means)

rownames(bird.means) <- NULL
colnames(bird.means) <- c("binomial","bioclim1", "bioclim2", "bioclim3", "bioclim4", "bioclim5", "bioclim6", "bioclim7", "bioclim8", "bioclim9", "bioclim10", "bioclim11", "bioclim12", "bioclim13", "bioclim14", "bioclim15", "bioclim16", "bioclim17", "bioclim18", "bioclim19")

#Write csv
write.csv(bird.means, 'data/bird.means.csv')

4) This process was repeated when we used gridded data for the LGM, LIG, NPP and human population density, while range size was obtained from the following code:

bird.range.size <- sapply(bird.ranges,
                       function(x) {(area(x))})
bird.range.size <- sapply(bird.range.size, function(x) {sum(x)})
bird.range.size <- as.data.frame(bird.range.size)
bird.range.size <- cbind.data.frame(shps.jetz, bird.range.size)
rownames(bird.range.size) <- NULL
colnames(bird.range.size) <- c("binomial", "range.size.m2")

#write.csv
write.csv(bird.range.size, 'data/bird.range.size.csv')

Generating biologically relevent predictors of extinction

Here we assess the relationship between sexual selection and extinction risk. However, in doing so we attempt to take into account as many other predictors of extinction as possible, primarily through environmental variables. Basically, our model structure seeks to contain:

Extinction/Diversification ~ Sexual selection
+ Range size
+ Short temporal variability of temperature (mean BIOCLIM4)
+ Spatial variability of temperature (PCA1) [residual.PCA1]
+ Long-term variability of temperature (LIG)
+ NPP

These predictors can be obtained from the compiled datasets seen here:

plumage.scores <- read.csv('data/plumage_scores.csv')
#Generate sexual dichromatism score: 
plumage.scores$SDi <- abs(plumage.scores$Male_plumage_score - plumage.scores$Female_plumage_score)

#Make DF with enviro variables and plumage scores
complete.dataframe <- read.csv('data/complete.dataframe.csv')
complete.dataframe <- left_join(plumage.scores %>% dplyr::select(binomial, Male_plumage_score, Female_plumage_score, SDi), complete.dataframe, by = "binomial")

Spatial Variability PCAs

To obtain estimates of spatial variability we can obtain PCAs of the standard errors of the bioclimatic variables (excluding seasonality in temperature and precipitation). From there we can run GAM models with the PCA components and range size and extract residuals. We do this as we expect variability will increase alongside increased range size, thus taking the residuals allows us to correct for this association.

Table S1: Loadings for the first two PCs, when we conduct a PCA on the variation (standard error, where n ~ 1,000) pf each bioclim variable except the standard error of temperature seasonality (se.bioclim4) and the standard error of precipitation seasonality (se.bioclim15), as these are already measures of variability and based on other bioclim variables. We use the standard error as, while we attempted to sample each species range 1,000 times, some points returned NA values and would thus affect the standard deviation, this inconvenience should not affect standard error.

restricted.data <- complete.dataframe %>% drop_na(bioclim1) #Drop species without environmental data
PCA.bioclim <- prcomp(restricted.data[c('se.bioclim1', 'se.bioclim2', 'se.bioclim3', 'se.bioclim5', 'se.bioclim6', 'se.bioclim7', 'se.bioclim8', 'se.bioclim9', 'se.bioclim10', 'se.bioclim11', 'se.bioclim12', 'se.bioclim13', 'se.bioclim14', 'se.bioclim16', 'se.bioclim17', 'se.bioclim18', 'se.bioclim19')], #Select se.bioclim
                      scale = TRUE, center = TRUE)
PCA.predictions <- predict(PCA.bioclim)
restricted.data <- cbind(restricted.data, PCA.predictions)

#Create table with PC1 and PC2
as.data.frame(PCA.bioclim$rotation[,1:2]) %>% `colnames<-`(c("PC1", "PC2")) %>% pander()

	PC1	PC2
se.bioclim1	-0.3328	0.08264
se.bioclim2	-0.2383	0.05458
se.bioclim3	-0.2317	-0.06584
se.bioclim5	-0.312	0.08091
se.bioclim6	-0.3288	0.08275
se.bioclim7	-0.2934	0.09748
se.bioclim8	-0.2916	0.1181
se.bioclim9	-0.3187	0.1125
se.bioclim10	-0.3111	0.06704
se.bioclim11	-0.3296	0.08775
se.bioclim12	-0.1386	-0.4246
se.bioclim13	-0.1596	-0.3083
se.bioclim14	-0.01133	-0.3857
se.bioclim16	-0.1708	-0.3254
se.bioclim17	-0.01923	-0.3932
se.bioclim18	-0.1433	-0.3066
se.bioclim19	-0.01797	-0.3815

#run GAM
PC1.gam <- mgcv::gam(PC1*(-1) ~ s(log(range.size.m2)), data = restricted.data, family = "gaussian") #Inverse as the PC1 loads to the negative (COUNTER-NTUITIVE)

#take residuals
restricted.data$residuals.PC1 <- residuals.gam(PC1.gam) 

#We can do the same for PC2
PC2.gam <- mgcv::gam(PC2 ~ s(log(range.size.m2)), data = restricted.data, family = "gaussian")

#Take residuals
restricted.data$residuals.PC2 <- residuals.gam(PC2.gam)

#plot
par(mfrow = c(1,2))
plot.gam(PC1.gam, residuals = T, main = 'PC1 residuls (temp)')
plot.gam(PC2.gam, residuals = T, main = 'PC2 residuals (precip)')

Figure S1: The relationship between spatial variability in the temperature components (PC1) and log-range size is relatively strong. But not as strong for spatial variability in precipitation (PC2). In the analysis we only used the residuals from PC1

Long term climate variability (LIG anomaly)

To gain estimates for change in climate over the past ~130,000 years we can use the difference in bioclim variables between the LIG and present values. The plots show the two PCs, broadly representing temperature and prescipitation. Data has also been provided for the last-glacial maximum (LGM) but was not used in the analysis.

Table S2: Loadings for the first two PCs of each of the PCAs for the difference in bioclimatic variables between today and the LIG. Here, PC1 is more heavily loaded for absoluted temperature difference, while PC2 is more heavily loaded for absolute difference in precipitation.

historical.variation.data <- as.data.frame(restricted.data[1])

#FOR LIG

historical.variation.data$bio1.LIG.diff <- abs(restricted.data$bioclim1 - restricted.data$LIG.bi1)
historical.variation.data$bio2.LIG.diff <- abs(restricted.data$bioclim2 - restricted.data$LIG.bi2)
historical.variation.data$bio3.LIG.diff <- abs(restricted.data$bioclim3 - restricted.data$LIG.bi3)
historical.variation.data$bio4.LIG.diff <- abs(restricted.data$bioclim4 - restricted.data$LIG.bi4)
historical.variation.data$bio5.LIG.diff <- abs(restricted.data$bioclim5 - restricted.data$LIG.bi5)
historical.variation.data$bio6.LIG.diff <- abs(restricted.data$bioclim6 - restricted.data$LIG.bi6)
historical.variation.data$bio7.LIG.diff <- abs(restricted.data$bioclim7 - restricted.data$LIG.bi7)
historical.variation.data$bio8.LIG.diff <- abs(restricted.data$bioclim8 - restricted.data$LIG.bi8)
historical.variation.data$bio9.LIG.diff <- abs(restricted.data$bioclim9 - restricted.data$LIG.bi9)
historical.variation.data$bio10.LIG.diff <- abs(restricted.data$bioclim10 - restricted.data$LIG.bi10)
historical.variation.data$bio11.LIG.diff <- abs(restricted.data$bioclim11 - restricted.data$LIG.bi11)
historical.variation.data$bio12.LIG.diff <- abs(restricted.data$bioclim12 - restricted.data$LIG.bi12)
historical.variation.data$bio13.LIG.diff <- abs(restricted.data$bioclim13 - restricted.data$LIG.bi13)
historical.variation.data$bio14.LIG.diff <- abs(restricted.data$bioclim14 - restricted.data$LIG.bi14)
historical.variation.data$bio15.LIG.diff <- abs(restricted.data$bioclim15 - restricted.data$LIG.bi15)
historical.variation.data$bio16.LIG.diff <- abs(restricted.data$bioclim16 - restricted.data$LIG.bi16)
historical.variation.data$bio17.LIG.diff <- abs(restricted.data$bioclim17 - restricted.data$LIG.bi17)
historical.variation.data$bio18.LIG.diff <- abs(restricted.data$bioclim18 - restricted.data$LIG.bi18)
historical.variation.data$bio19.LIG.diff <- abs(restricted.data$bioclim19 - restricted.data$LIG.bi19)

historical.variation.data <- historical.variation.data %>% drop_na(bio1.LIG.diff)

#Run a PCA of the difference removing the variation variables (4 and 15)

#For LIG
PCA.LIG.bioclim <- prcomp(historical.variation.data[c(
'bio1.LIG.diff',
'bio2.LIG.diff',
'bio3.LIG.diff',
'bio5.LIG.diff',
'bio6.LIG.diff',
'bio7.LIG.diff',
'bio8.LIG.diff',
'bio9.LIG.diff',
'bio10.LIG.diff',
'bio11.LIG.diff',
'bio12.LIG.diff',
'bio13.LIG.diff',
'bio14.LIG.diff',
'bio16.LIG.diff',
'bio17.LIG.diff',
'bio18.LIG.diff',
'bio19.LIG.diff')],
                      scale = TRUE, center = TRUE)

#Create table with PC1 and PC2
as.data.frame(PCA.LIG.bioclim$rotation[,1:2]) %>% `colnames<-`(c("PC1.LIG", "PC2.LIG")) %>% pander()

	PC1.LIG	PC2.LIG
bio1.LIG.diff	-0.2956	-0.01943
bio2.LIG.diff	-0.2442	0.0745
bio3.LIG.diff	-0.1436	0.2085
bio5.LIG.diff	-0.2867	-0.2566
bio6.LIG.diff	-0.4094	0.03346
bio7.LIG.diff	-0.4012	-0.09706
bio8.LIG.diff	0.07108	-0.2385
bio9.LIG.diff	-0.3365	-0.00753
bio10.LIG.diff	-0.07014	-0.3766
bio11.LIG.diff	-0.4201	-0.004
bio12.LIG.diff	-0.05181	0.3337
bio13.LIG.diff	-0.2552	0.2591
bio14.LIG.diff	0.05701	0.3298
bio16.LIG.diff	-0.1946	0.3164
bio17.LIG.diff	0.03896	0.3603
bio18.LIG.diff	0.06173	0.277
bio19.LIG.diff	0.08743	0.2856

#Now we can predict the PCA results: 
PCA.LIG.predictions <- as.data.frame(predict(PCA.LIG.bioclim)*(-1)) #So that higher numbers mean more variation
PCA.LIG.predictions <- rename(PCA.LIG.predictions, PC1.LIG = PC1,
       PC2.LIG = PC2)

#Bind them to dataframe, taking only the first two PCAs
historical.variation.data <- cbind(historical.variation.data, PCA.LIG.predictions[1:2])

#Now to the restricted dataframe
restricted.data <- right_join(restricted.data, historical.variation.data %>% dplyr::select(binomial, PC1.LIG, PC2.LIG), by = 'binomial')

PCA.plot <- historical.variation.data %>% ggplot(aes(x = PC1.LIG, y= PC2.LIG))+
  geom_point(shape = 21)+
  theme_minimal()

PCA.plot.m <- ggExtra::ggMarginal(
  p = PCA.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = 'gray'
)

grid.arrange(PCA.plot.m)

Figure S2: The distribution of the two PCA compontents for the absolute difference in bioclimatic variables between today and the LIG.

Correlations between environmental predictors

We can check the correlations between environmental predictors that we will use in our model. Specifically we want to know whether the following environmental variables are correlated:

Range size
Short temporal variability of temperature (mean BIOCLIM4)
Spatial variability of temperature (PCA1) [residual.PCA1]
Long-term variability of temperature (LIG)
NPP

To check some of the correlations between the environmental predictors to be used in the PGLS models we can inspect the following correlation plot with the correlation value plotted:

#Draw plot
restricted.data$log.range.size <- log(restricted.data$range.size.m2)
pairs(restricted.data[,c("log.range.size", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP")], lower.panel=panel.smooth, upper.panel=panel.cor)

Figure S3: The highest correlation is between Long term temperature variation and seasonal variation. However this correlation is not excessively high (0.53) so we can be satisfied in including these predictors within our model.

Sexual Dimorphism

The sexual dimorphism dataset looks to have an overdispersed distribution (see phylogenic plot in manuscript). Unfortunately, transformations do not greatly improve the distribution. Our model fit is expected to be reduced by this but further transformations are not obvious and would risk problems of interperetation.

Analysis

Simple speciation measures: Diversification rate (DR) and Node Density (ND)

Diversification rates can be simply estimated from thew ‘tippiness of the phylogeny’. Here a log equal splits (also called DR) value can be achieved that is weighted to more recent diversification. It was used by Jetz et al. (2012) and recently (in comparison with BAMM) in a study ivestigation bird diversity alongside elevation (Quintero and Jetz 2018). We can look at speciation rate using ES-Sim adapted by Harvey Michael et al. (2017) and found here https://github.com/mgharvey/ES-sim.

The functions that generate these tip rates estimates are: calculate.log.es, calculate.nd and calculate.tb. Here they are run on 100 trees plus a full MCC tree of passerine birds based on 2500 samples of the posterior (Jetz et al. 2012). This was done using the phanghorn (Schliep 2010).

#Obtain DR (log(es)) estimates, by calling the calculate.dr function
passerine.trees <- read.nexus('bamm_extinction/data/trees/Passerine_Trees_Full.nex')

#DR, ND and TB can be calculated on the 100 trees by: 
  
# es.list <- sapply(passerine.trees, calculate.log.es)
# nd.list <- sapply(passerine.trees, calculate.nd)
# tb.list <- sapply(passerine.trees, calculate.tb)

#To avoid re-running the time-consuming rates, we can save them and load them: 
es.list <- readRDS("data/es.list.rds")
nd.list <- readRDS("data/nd.list.rds")
tb.list <- readRDS("data/tb.list.rds")

#Also load in our MCC tree
MCC.passerine <- read.tree('data/MCC.passerine.tre') #Tree based on 2500 samples of posterior

To gain an estimate of variation between trees we can take the values and obtain a mean, variance and coeeficient of variation (CV). We obtain a CV value to account for unequal variation at higher estimates. The CV is expressed as the ratio of the standard deviation (\(\sigma\)) to the mean (\(\mu\)) (converting to a percentage not needed):

\(C_V = \frac{\sigma}{\mu} \times 100\)

Furthermore, we can account for overdispersion by obtaining a transformed version of CV. In such a case we can calculate CV for a log-normal distribution:

\({\widehat {c_{\rm {v}}}}_{\rm {ln}}={\sqrt {\mathrm {e} ^{s_{\rm {ln}}^{2}}-1}}\)

Given that our ES value is log-transformed, the calculation of \({\widehat {c_{\rm {v}}}}_{\rm {ln}}\) is perhaps more ideal. As for ND we can see from the a qqnorm plot that value of ND is normally distributed, therefore the regular CV is more appropriate for this distribution.

#For ND
#transpose the list so that the elements are the species
t.nd.list <- nd.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
nd.values <- as.data.frame(t(sapply(t.nd.list, function(x) sapply(x, function(x) (x)))))
nd.values <- nd.values %>% tibble::rownames_to_column("binomial")
nd.summary <- melt(nd.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.nd = mean(value), #use normal mean, not a rate
  var.nd = var(value), #Don't really need this if we use CV
  CV.nd = (sd(value)/mean(value)) #Given normal distribution of ND, standard CV formula is appropriate
)

#For ES
#transpose the list so that the elements are the species
t.es.list <- es.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
es.values <- as.data.frame(t(sapply(t.es.list, function(x) sapply(x, function(x) (x)))))
es.values <- es.values %>% tibble::rownames_to_column("binomial")
es.summary <- melt(es.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.loges = log(1/mean(1/exp(value))), #use harmonic mean for rates and log-transform post calculation of HM
  var.loges = var(value), #Don't really need this if we use CV
  CV.loges = sqrt(exp(var(value)) - 1) #CV for log-normal distribution of ES
)

#For TB
#transpose the list so that the elements are the species
t.tb.list <- tb.list %>% purrr::transpose()
#not too sure why this works but turn the nested list into a dataframe
tb.values <- as.data.frame(t(sapply(t.tb.list, function(x) sapply(x, function(x) (x)))))
tb.values <- tb.values %>% tibble::rownames_to_column("binomial")
tb.summary <- melt(tb.values, id.vars = "binomial") %>% group_by(binomial) %>% summarise(
  mean.tb = mean(value), #use normal mean; not a rate
  var.tb = var(value), #Don't really need this if we use CV
  CV.tb = (sd(value)/mean(value)) #Given normal distribution of tb, standard CV formula is appropriate
)

MCC.nd.df <- as.data.frame(calculate.nd(MCC.passerine)) %>% tibble::rownames_to_column("binomial")
MCC.dr.df <- as.data.frame(calculate.log.es(MCC.passerine)) %>% tibble::rownames_to_column("binomial")

#Join the dataframes
dr.summary <- full_join(nd.summary, es.summary, by = "binomial")
dr.summary <- full_join(dr.summary, tb.summary, by = "binomial")
dr.summary <- full_join(dr.summary, MCC.nd.df, by = "binomial")
dr.summary <- full_join(dr.summary, MCC.dr.df, by = "binomial")

dr.summary <- dr.summary %>% rename(TipLabel = binomial, MCC.ND = nd, MCC.DR = es)

saveRDS(dr.summary, 'data/dr.summary.rds')

dr.summary <- readRDS('data/dr.summary.rds')
#Plot a summary of logES and TB with the DRs with their weightings
dr.plot <- dr.summary %>% ggplot(aes(x = mean.loges, y = mean.nd, 
                          fill = (1/(CV.nd))/sum((1/(CV.nd)))*100, 
                      size = (1/(CV.loges))/sum((1/(CV.loges)))*100))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = mean.loges - 0.674*sqrt(var.loges), xmax = mean.loges + 0.674*sqrt(var.loges)), 
                 size = 0.015)+
  geom_errorbar(aes(ymin = mean.nd - 0.674*sqrt(var.nd), ymax = mean.nd + 0.674*sqrt(var.nd)), 
                size = 0.015)+
  theme(legend.position="bottom")+
  labs(size = 'Weight (%) [using DR]', y='Node Density [ND]', x= 'Log Diversification Rate [DR]', fill = 'Weight (%) [using ND] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "RdPu", direction = 1, breaks=c(0,0.015,.03),
                           limits=c(0,0.03))
marginal.dr <-ggMarginal(
  p = dr.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = '#0000009C',
  fill = '#D12E6769'
)

grid.newpage()
grid.draw(marginal.dr)

Figure S4: \(\lambda_{DR}\) and \(\lambda_{ND}\) are correlated in their estimation of species rate. Furthermore, we see similar variation across each species in the estimates across 100 trees as the inverse of the coefficeint of variation (weight) is approximately equivalent across all values of ND and DR (colours and size). Notably, these estimates speciation rate are relatively uncertain across all species, with 50 % CIs in both axes plotted.

Combine the estimates of DR and ND with the restricted dataframe containing all the other variables. We can inspect the trend in SDi in relation to the speciation rates.

restricted.data <- left_join(restricted.data, dr.summary, by = 'TipLabel')

es.plot <- restricted.data %>% ggplot(aes(y = mean.loges, x = SDi, fill = mean.loges))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "YlOrBr", direction = 1, guide = FALSE)+
ylab("mean.DR")+
theme_minimal()

nd.plot <- restricted.data %>% ggplot(aes(y = mean.nd, x = SDi, fill = mean.nd))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "Greens", direction = 1, guide = FALSE)+
theme_minimal()

tb.plot <- restricted.data %>% ggplot(aes(y = mean.tb, x = SDi, fill = mean.tb))+
geom_point(shape = 21, size = 1.5)+
geom_smooth(show.legend = FALSE, color = "grey20", method = "lm")+
scale_fill_distiller(palette = "PuBu", direction = 1, guide = FALSE)+
theme_minimal()

grid.arrange(es.plot, nd.plot, tb.plot, nrow = 3)

Figure S5: Scatter plots showing the raw relationship between sexual dichromatism (SDi) and speciation rates. Across all three measures the pattern and spread is similar, with no obvious relationship. Note that terminal branch lenght (mean.tb) was not used in the analysis.

BAMM measures of speciation and extinction

Set up BAMM parameters

For the use of BAMM we used the following code to generate the parameters across the 100 trees. Each paramaeter value is specified in this code chunk to ensure reproducibility when re-running on any of the trees. These same parameters were also used for the MCC tree (with the seed set at 2500).

name.passerine.tree <- names(passerine.trees)

priors <- sapply(name.passerine.tree, function(x) {
  setBAMMpriors(passerine.trees[[x]], outfile = NULL)
})

sapply(name.passerine.tree, function(x) {
  write.tree(passerine.trees[[x]], paste("data/bamm_files/", x, ".tre", sep=""))
})

# Here is a block of parameters for the control file. We can make a control file for each tree:
params <- list()
for (x in name.passerine.tree) {

# GENERAL SETUP AND DATA INPUT

params[[x]] <- list(modeltype = 'speciationextinction',
# Specify "speciationextinction" or "trait" analysis
                                  
treefile = paste(x, ".tre", sep=""),
# File name of the phylogenetic tree to be analyzed

runInfoFilename = 'run_info.txt',
# File name to output general information about this run

sampleFromPriorOnly = 0,
# Whether to perform analysis sampling from prior only (no likelihoods computed)

runMCMC = 1,
# Whether to perform the MCMC simulation. If runMCMC = 0, the program will only
# check whether the data file can be read and the initial likelihood computed

loadEventData = 0,                       
# Whether to load a previous event data file

eventDataInfile = 'event_data_in.txt',
# File name of the event data file to load, used only if loadEventData = 1

initializeModel = 1,
# Whether to initialize (but not run) the MCMC. If initializeModel = 0, the
# program will only ensure that the data files (e.g., treefile) can be read

useGlobalSamplingProbability = 1,
# Whether to use a "global" sampling probability. If False (0), expects a file
# name for species-specific sampling probabilities (see sampleProbsFilename)
                                        
globalSamplingFraction = 1,
# The sampling probability. If useGlobalSamplingProbability = 0, this is ignored
# and BAMM looks for a file name with species-specific sampling fractions

sampleProbsFilename = 'sample_probs.txt',
# File name containing species-specific sampling fractions

seed = as.numeric(gsub("tree_", "", x, perl = TRUE)),
# Seed for the random number generator. Set for reproducibility to the number of the treefile

overwrite = 1,
# If True (1), the program will overwrite any output files in the current
# directory (if present)


# PRIORS

expectedNumberOfShifts = 100,
# prior on the number of shifts in diversification
# Suggested values: 
#     expectedNumberOfShifts = 1.0 for small trees (< 500 tips)
#  expectedNumberOfShifts = 10 or even 50 for large trees (> 5000 tips) 
 
lambdaInitPrior = as.numeric(priors['lambdaInitPrior', x]),
# Prior (rate parameter of exponential) on the initial lambda value for rate
# regimes

lambdaShiftPrior = 0.05,
# Prior (std dev of normal) on lambda shift parameter for rate regimes
# You cannot adjust the mean of this distribution (fixed at zero, which is
# equal to a constant rate diversification process)

muInitPrior = as.numeric(priors['muInitPrior', x]),
# Prior (rate parameter of exponential) on extinction rates  

lambdaIsTimeVariablePrior = 1,
# Prior (probability) of the time mode being time-variable (vs. time-constant)
            

# MCMC SIMULATION SETTINGS & OUTPUT OPTIONS

numberOfGenerations = '100000000',
# Number of generations to perform MCMC simulation

mcmcOutfile = 'mcmc_out.txt',
# File name for the MCMC output, which only includes summary information about
# MCMC simulation (e.g., log-likelihoods, log-prior, number of processes)

mcmcWriteFreq = 1000,
# Frequency in which to write the MCMC output to a file

eventDataOutfile = 'event_data.txt',
# The raw event data (these are the main results). ALL of the results are
# contained in this file, and all branch-specific speciation rates, shift
# positions, marginal distributions etc can be reconstructed from this output.
# See R package BAMMtools for working with this output

eventDataWriteFreq = 1000,
# Frequency in which to write the event data to a file

printFreq = 10000,
# Frequency in which to print MCMC status to the screen

acceptanceResetFreq = 1000,
# Frequency in which to reset the acceptance rate calculation
# The acceptance rate is output to both the MCMC data file and the screen

outName = x,
# Optional name that will be prefixed on all output files (separated with "_")
# If commented out, no prefix will be used


# OPERATORS: MCMC SCALING OPERATORS

updateLambdaInitScale = 2,
# Scale parameter for updating the initial speciation rate for each process

updateLambdaShiftScale = 0.1,
# Scale parameter for the exponential change parameter for speciation

updateMuInitScale = 2,
# Scale parameter for updating initial extinction rate for each process

updateEventLocationScale = 0.1,
# Scale parameter for updating LOCAL moves of events on the tree
# This defines the width of the sliding window proposal
 
updateEventRateScale = 4,
# Scale parameter (proportional shrinking/expanding) for updating
# the rate parameter of the Poisson process

# OPERATORS: MCMC MOVE FREQUENCIES

updateRateEventNumber = 1,
# Relative frequency of MCMC moves that change the number of events

updateRateEventPosition = 0.25,
# Relative frequency of MCMC moves that change the location of an event on the
# tree

updateRateEventRate = 1,
# Relative frequency of MCMC moves that change the rate at which events occur 

updateRateLambda0 = 1,
# Relative frequency of MCMC moves that change the initial speciation rate
# associated with an event

updateRateLambdaShift = 1,
# Relative frequency of MCMC moves that change the exponential shift parameter
# of the speciation rate associated with an event

updateRateMu0 = 1,
# Relative frequency of MCMC moves that change the extinction rate for a given
# event

updateRateLambdaTimeMode = 0,
# Relative frequency of MCMC moves that flip the time mode
# (time-constant <=> time-variable)

localGlobalMoveRatio = 10,
# Ratio of local to global moves of events 


# INITIAL PARAMETER VALUES

lambdaInit0 = 0.032,
# Initial speciation rate (at the root of the tree)

lambdaShift0 = 0,
# Initial shift parameter for the root process

muInit0 = 0.005,
# Initial value of extinction (at the root)

initialNumberEvents = 0,
# Initial number of non-root processes


# METROPOLIS COUPLED MCMC

numberOfChains = 1,
# Number of Markov chains to run

deltaT = 0.01,
# Temperature increment parameter. This value should be > 0
# The temperature for the i-th chain is computed as 1 / [1 + deltaT * (i - 1)]

swapPeriod = 1000,
# Number of generations in which to propose a chain swap

chainSwapFileName = 'chain_swap.txt',
# File name in which to output data about each chain swap proposal.
# The format of each line is [generation],[rank_1],[rank_2],[swap_accepted]
# where [generation] is the generation in which the swap proposal was made,
# [rank_1] and [rank_2] are the chains that were chosen, and [swap_accepted] is
# whether the swap was made. The cold chain has a rank of 1.


# NUMERICAL AND OTHER PARAMETERS

minCladeSizeForShift = 3,
# Allows you to constrain location of possible rate-change events to occur
# only on branches with at least this many descendant tips. A value of 1
# allows shifts to occur on all branches. 

segLength = 0.025,
# Controls the "grain" of the likelihood calculations. Approximates the
# continuous-time change in diversification rates by breaking each branch into
# a constant-rate diversification segments, with each segment given a length
# determined by segLength. segLength is in units of the root-to-tip distance of
# the tree. So, if the segLength parameter is 0.01, and the crown age of your
# tree is 50, the "step size" of the constant rate approximation will be 0.5.
# If the value is greater than the branch length (e.g., you have a branch of
# length < 0.5 in the preceding example) BAMM will not break the branch into
# segments but use the mean rate across the entire branch.

outName = x)
  }

bammcontrolfile <- list()
for (x in name.passerine.tree) {
  bammcontrolfile[x] <- paste("data/bamm_files/control_", x, ".txt", sep="")
}

# Now writing control parameters to file

for (x in name.passerine.tree) {generateControlFile(file = bammcontrolfile[[x]], type = "diversification", params = params[[x]])}

Run analysis

BAMM can be run through the terminal through the following syntax: bamm -c control_tree_xxxx.txt. To generate these commands we can use a loop function, from which we get:

bamm.commands <- list()
for (x in name.passerine.tree) {
  bamm.commands[x] <- paste("bamm -c control_", x, ".txt", sep="")
}

The analysis was run over multiple CPU’s, each generating a respective MCMC and EventData output. Due to the size of the event data file (~ 50 Gb in total) they are not included as supplementary material here. However we have simplified the event data objects into tip rate estimates of the mean and variance across 100 trees + MCC.

Read in the event data and extract tip data

Now we have a series of trees and event data we can read in. Firstly and most importantly let’s check for convergence in the MCC tree through the following:

Table S3: Effective sample size (ESS) for the two key BAMM parameters (number of evolutionary shifts and log-Likelihood) for the run on the MCC indicate that BAMM converges (ESS > 200).

#Read in the tree and MCMC to check for convergence 
MCC.BAMM.tree  <- read.tree("data/BAMM_MCC/MCC.passerine.tre") #Same as other MCC tree already loaded
MCC.BAMM.mcmc <- read.csv( "data/BAMM_MCC/tree_MCC_mcmc_out.txt" , stringsAsFactors=F)

#Plot of convergence can be generated by:
#plot(MCC.BAMM.mcmc$logLik ~ MCC.BAMM.mcmc$generation)

#Looks like it has converged so let's discard burn in: 
burnstart <- floor(0.1 * nrow(MCC.BAMM.mcmc))
postburn <- MCC.BAMM.mcmc[burnstart:nrow(MCC.BAMM.mcmc), ]

#We can also check effective population sizes of the log-likelihhod and number of shift events in each sample
#We want at least 200 (although that's on the low side)

cbind(effectiveSize(postburn$N_shifts), effectiveSize(postburn$logLik)) %>% `colnames<-`(c("N_Shifts", "logLik")) %>% `rownames<-`("Effective Sample Size") %>% pander()

	N_Shifts	logLik
Effective Sample Size	2456	1302

Those values are well above 200, so we are good to use the MCC data in analysis! We can also check the convergence of BAMM across 100 runs of BAMM. The Raw MCMCs are not included in this file or repository but we can read in a dataframe that has extracted the ESS for each of the runs.

Table S4: For the 100 trees that BAMM was run on effective sample size (ESS) for the two key BAMM parameters (number of evolutionary shifts and log-Likelihood) also indicates that BAMM converges with the minimum for each parameter across the 100 trees being over 200.

ESS <- readRDS('data/ESS.rds')
summary(ESS) %>% pander()

N_Shifts	logLik
Min. : 723.4	Min. : 278.7
1st Qu.:1396.1	1st Qu.: 511.3
Median :1618.6	Median : 638.6
Mean :1709.5	Mean : 669.9
3rd Qu.:2070.8	3rd Qu.: 804.8
Max. :3425.6	Max. :1244.1

Given that it appears BAMM converges we can we can make a data frame with the mean and variance for extinction and speciation tip rates from the large event data set for the MCC with the following code and then plot the variation we see in tip-rates.

# Read in Event Data
MCC.BAMM.ED <- getEventData(MCC.BAMM.tree,  "data/BAMM_MCC/tree_MCC_event_data.txt", burnin=0.1, nsamples=1000)

## Reading event datafile:  data/BAMM_MCC/tree_MCC_event_data.txt 
##      ...........
## Read a total of 100000 samples from posterior
## 
## Discarded as burnin: GENERATIONS <  9999000
## Analyzing  1000  samples from posterior
## 
## Setting recursive sequence on tree...
## 
## Done with recursive sequence

saveRDS(MCC.BAMM.ED, 'data/MCC.BAMM.ED.rds')
#From the Event Data we can extract 

library(purrr)
#BAMM.EventData <- readRDS('data/BAMM.EventData.rds')

#Big lapply over each tree in BAMM event data
BAMM.extraction.function <- function(x) {
######Get mean and var for lambda
#Transpose list so each element in the list is a species
transposed.lambda <- lapply(purrr::transpose(x$tipLambda), unlist)

#Now turn it into a df with mean and variance
lambda <- sapply(transposed.lambda, function(x) {
  mean.lambda = mean(x)
  var.lambda = var(x)
  return(c(mean.lambda, var.lambda))
}) %>% t() %>% as.data.frame() %>% `colnames<-`(c("mean.lambda", "var.lambda"))

lambda$TipLabel <- x[["tip.label"]]

#####NOW FOR Extinction

#Transpose list so each element in the list is a species
transposed.mu <- lapply(purrr::transpose(x$tipMu), unlist)

#Now turn it into a df with mean and variance
mu <- sapply(transposed.mu, function(x) {
  mean.mu = mean(x)
  var.mu = var(x)
  return(c(mean.mu, var.mu))
}) %>% t() %>% as.data.frame() %>% `colnames<-`(c("mean.mu", "var.mu"))

mu$TipLabel <- x[["tip.label"]]

left_join(lambda, mu, by = "TipLabel")
}

MCC.BAMM.df <- BAMM.extraction.function(MCC.BAMM.ED)

#Save df for later use
saveRDS(MCC.BAMM.df, 'data/MCC.BAMM.df.rds')

MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')
#Plot a summary of logES and TB with the DRs with their weightings
BAMM.MCC.plot <- MCC.BAMM.df %>% ggplot(aes(x = log(mean.lambda), y = log(mean.mu), 
                          fill = (1/(log(var.lambda)))/sum((1/(log(var.lambda))))*100, 
                      size = (1/(log(var.mu)))/sum((1/(log(var.mu))))*100))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = log(mean.lambda) - 0.674*log(sqrt(var.lambda)), xmax = log(mean.lambda) + 0.674*log(sqrt(var.lambda))), 
                 size = 0.0025)+
  geom_errorbar(aes(ymin = log(mean.mu) - 0.674*log(sqrt(var.mu)), ymax = log(mean.mu) + 0.674*log(sqrt(var.mu))), 
                size = 0.0025)+
  xlim(-8, 3)+
  ylim(-10,2)+
  # scale_x_continuous(trans = "log")+
  # scale_y_continuous(trans = "log")+
  theme(legend.position="bottom")+
  labs(size = 'Weight (%) [using Lambda]', y='Log Extinction [Mu]', x= 'Log Speciation [Lambda]', fill = 'Weight (%) [using Mu] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "Reds", direction = 1, breaks=c(0,.15,0.3),
                           limits=c(-0,0.3))

BAMM.variance <- ggExtra::ggMarginal(
  p = BAMM.MCC.plot,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = '#BA3B1C91'
)

grid.newpage()
grid.draw(BAMM.variance)

Figure S5: The tip-rate estimate for BAMM are highly variable within each run of BAMM. Across most species there is high variability in the posterior distribution of tip-rate estimates. Here we show mean values, weights based on the variance and 50 % CIs.

Analysis of BAMM results

Given the high variability in tip-rate estimates from the above plot, below we peform some diagnostics on BAMM to demonstrate that the variability is unlikely an error in sampling or parameters, rather an inherent aspect of BAMM. The following is not inherently necessary to understandfor the conclusions drawn from the paper, however they do raise a set of methodological questions about BAMM that warrant further investigation.

Credible number of shifts

To plot the credible shift set, we need the prior distribution on the number of rate shifts (this is generated internally by BAMMtools). We can then estimate the credible set of rate shifts using the BAMMtools function credibleShiftSet:

css <- credibleShiftSet(MCC.BAMM.ED, expectedNumberOfShifts=100, threshold=5, set.limit = 0.95)

#Now we obtain the number of distinct shifts: (Out of 1000 samples this is super high, essentially each one is distinct)

css$number.distinct

## [1] 950

This value indicates that from a posterior sample of 1000 we have 950 unique rate shifts, which is very high. Arguably, across such a large phylogenetic tree, even if the number of shifts reaches convergence the location of those > 50 shifts will be quite variable, with the variation likely impacting tip-rate estimates. Even if we increase the sampling of the posterior to 2,500 instead of 1,000 we still see that the CSS is high (close to 2500). We can have a look what our top few shifts are:

summary(css)

## 
##  95 % credible set of rate shift configurations sampled with BAMM
## 
## Distinct shift configurations in credible set:  950
## 
## Frequency of 9 shift configurations with highest posterior probability:
## 
## 
##    rank     probability cumulative  Core_shifts
##          1      0.001      0.001         54
##          2      0.001      0.002         44
##          3      0.001      0.003         46
##          4      0.001      0.004         47
##          5      0.001      0.005         49
##          6      0.001      0.006         45
##          7      0.001      0.007         47
##          8      0.001      0.008         43
##          9      0.001      0.009         45
## 
## ...omitted 941 additional distinct shift configurations
## from the credible set. You can access the full set from your 
## credibleshiftset object

Notably, each of the top shifts has low probability, indicating that out of the entire posterior sample we cannot differentiate which shift configuration is more likely. Based on the follwing eastimates of the BayesFactor, we are convinced that the number of shifts is non-zero, however it is still quite a wide distribution, sourced from potentially weakly non-informative priors. But, this is at odds with the ESS and the following plot:

round(computeBayesFactors(MCC.BAMM.mcmc, expectedNumberOfShifts=100, burnin=0.1)[,1], digits = 2)

##      42      43      44      45      46      47      48      49      50 
##    1.00   17.17   24.48   49.45  120.71  283.77  454.33  784.80 1291.85 
##      51      52      53      54      55      56      57      58      59 
## 1813.33 2561.62 3608.07 4574.91 5492.44 6322.11 7141.12 7720.26 7886.28 
##      60      61      62      63      64      65      66      67      68 
## 7846.73 7661.83 7116.15 6466.36 5749.34 4811.16 4136.80 3338.17 2674.70 
##      69      70      71      72      73      74      75      76      77 
## 2137.61 1593.48 1218.40  868.01  643.91  452.36  311.07  199.17  135.99 
##      78      79      80      81      82      83      84      85      87 
##   85.85   65.03   40.87   14.74   13.40    9.02   10.63    1.53    1.56

plotPrior(MCC.BAMM.mcmc, expectedNumberOfShifts=100)

Figure S6: The apparant convergence of the number of shifts in the posterior is at odds with the variability seen in the CSS and although there is greater certainty in the number of shifts being within the range above, the position of those shifts remain highly variable.

We can compare our run of BAMM against Harvey et al. (2017) who used BAMM on a genetic-only MCC tree with different parameters

load('data/Hackett_split_eventsample.rda')
css3 <- credibleShiftSet(ed, expectedNumberOfShifts=100, threshold=5, set.limit = 0.95)
css3$number.distinct

## [1] 1188

Based on Their posterior sample of 1250, 1188 are unique. Like our run on the MCC, we find high uncertainity in the position of shift-configurations.

Let’s run a model based on genetic data from the Harvey MCC

Harvey.BAMM.df <- BAMM.extraction.function(ed)

Harvey.BAMM.df %>% ggplot(aes(x = log(mean.lambda), y = log(mean.mu), 
                          fill = (1/(log(var.lambda)))/sum((1/(log(var.lambda))))*100, 
                      size = (1/(log(var.mu)))/sum((1/(log(var.mu))))*100))+
  geom_point(shape = 21, alpha = 0.5)+
  theme_minimal()+
  geom_errorbarh(aes(xmin = log(mean.lambda) - 0.674*log(sqrt(var.lambda)), xmax = log(mean.lambda) + 0.674*log(sqrt(var.lambda))), 
                 size = 0.0025)+
  geom_errorbar(aes(ymin = log(mean.mu) - 0.674*log(sqrt(var.mu)), ymax = log(mean.mu) + 0.674*log(sqrt(var.mu))), 
                size = 0.0025)+
  xlim(-8, 3)+
  ylim(-12.5,1)+
  # scale_x_continuous(trans = "log")+
  # scale_y_continuous(trans = "log")+
  theme(legend.position="bottom")+
  labs(size = 'Weight (%) [using Lambda]', y='Log Extinction [Mu]', x= 'Log Speciation [Lambda]', fill = 'Weight (%) [using Mu] \n')+
  scale_fill_distiller(guide = "colorbar", palette = "Reds", direction = 1, breaks=c(0,.15,0.3),
                           limits=c(-0,0.3))

# BAMM.variance <- ggExtra::ggMarginal(
#   p = BAMM.MCC.plot,
#   type = 'density',
#   margins = 'both',
#   size = 5,
#   colour = 'black',
#   fill = '#BA3B1C91'
# )

Figure S7: Based on the event data from the BAMM run by Harvey et al. (2017) we find that in both our case and theirs the tip-rate estimates for many species are extrememly variable. here we present mean estimates with 50 % CIs.

PGLS Models

The method behind running the PGLS models is as follows:

1: Estimate the phylogenetic signal (\(\lambda\)) by running a model without interactions and all six predictor variables. The value of \(\lambda\) obtained here will be fixed in all successive models. The value is fixed for subsequent models as independently estimating it in each case becomes computationally intensive. Based on preliminary analysis even slight differences in lambda lead to large reductions in AIC.

2: Create a global model of six predictor variables plus the fiver interactions between sexual dichromatism and the other variables.

3: Dredge the global model but fix the six independent predictors, hence conducting model selection on the interaction terms.

4: Take the top model in the MCC model and run it on the 100 phylogenetic trees.

5: Repeat this step for DR, ND, BAMM-speciation and BAMM-extinction.

PGLS Models on DR and ND

Using corPagel we can estimate the phylogenetic signal for a model with all predictors (interactions do not appear to affect the estimate of \(\lambda\)):

#Prune tree
pruned.MCC.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(restricted.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(restricted.data) <- restricted.data$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.corPagel <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.tree, fixed = FALSE), 
                data = restricted.data, 
                method = "REML")
saveRDS(MCC.DR.corPagel, 'data/MCC.DR.corPagel.rds')

#Run a corPagel model to estimate lambda for ND
MCC.ND.corPagel <- gls(MCC.ND ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.tree, fixed = FALSE), 
                data = restricted.data, 
                method = "REML")
saveRDS(MCC.ND.corPagel, 'data/MCC.ND.corPagel.rds')

Inspect the \(\lambda\) value, which can then be fixed for successive models.

MCC.DR.corPagel <- readRDS('data/MCC.DR.corPagel.rds')
MCC.ND.corPagel <- readRDS('data/MCC.ND.corPagel.rds')
MCC.DR.corPagel[["modelStruct"]][["corStruct"]] %>% `names<-`("DR lambda") %>% pander()

DR lambda: 0.985

MCC.ND.corPagel[["modelStruct"]][["corStruct"]] %>% `names<-`("ND lambda") %>% pander()

ND lambda: 0.9996

The lambda is high and not too dissimilar than if we assume Brownian Motion \(\lambda = 1\). However, given the large sample size this difference may have an effect on the results so, to play safe we included it as a fixed value for \(\lambda\) in all successive models.

#Run model for DR
MCC.DR <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                correlation = corPagel(0.985, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

#Run model for ND
MCC.ND <- gls(MCC.ND ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP
                      + SDi*log(range.size.m2)
                      + SDi*bioclim4
                      + SDi*residuals.PC1
                      + SDi*PC1.LIG
                      + SDi*NPP,
                correlation = corPagel(0.9996, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("restricted.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

dredged.ND.model <- pdredge(MCC.ND, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust, trace = TRUE)

saveRDS(dredged.ND.model, "data/dredged.ND.model.rds")

dredged.DR.model <- pdredge(MCC.DR, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust, trace = TRUE)

saveRDS(dredged.DR.model, "data/dredged.DR.model.rds")

Table S5: The dredged models both show the top model is one with no interactions, with \(\delta AICc > 4\) in both cases. We can be resonably confident that interactions are unlikely to affect the pattern of speciation we see in passerine birds.

dredged.DR.model <- readRDS("data/dredged.DR.model.rds")
dredged.ND.model <- readRDS("data/dredged.ND.model.rds")
kable(dredged.DR.model, "html", caption = "DR Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

DR Dredge Table
	(Intercept)	bioclim4	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	SDi	bioclim4:SDi	log(range.size.m2):SDi	NPP:SDi	PC1.LIG:SDi	residuals.PC1:SDi	df	logLik	AICc	delta	weight
0	-2.879811	2.14e-05	-0.0065780	3e-07	0.0015864	-0.0024295	-0.0012787	NA	NA	NA	NA	NA	8	-3409.547	6835.118	0.00000	9.988652e-01
8	-2.879495	2.14e-05	-0.0065788	3e-07	-0.0001028	-0.0023765	-0.0012942	NA	NA	NA	0.0003487	NA	9	-3416.089	6850.208	15.09006	5.281325e-04
16	-2.879345	2.15e-05	-0.0066034	3e-07	0.0015509	-0.0015982	-0.0013082	NA	NA	NA	NA	-0.0001638	9	-3416.488	6851.007	15.88895	3.542138e-04
2	-2.884345	2.14e-05	-0.0064052	3e-07	0.0015991	-0.0024487	-0.0003035	NA	-0.0000369	NA	NA	NA	9	-3416.842	6851.714	16.59633	2.486904e-04
1	-2.878432	1.57e-05	-0.0065688	3e-07	0.0015872	-0.0024180	-0.0015875	1.2e-06	NA	NA	NA	NA	9	-3421.210	6860.452	25.33361	3.150533e-06
24	-2.878703	2.16e-05	-0.0066190	3e-07	-0.0004658	-0.0010549	-0.0013437	NA	NA	NA	0.0004121	-0.0002586	10	-3422.848	6865.734	30.61631	2.245227e-07
4	-2.884812	2.14e-05	-0.0065580	7e-07	0.0015560	-0.0024009	-0.0003587	NA	NA	-1e-07	NA	NA	9	-3424.307	6866.644	31.52618	1.424567e-07
10	-2.890455	2.13e-05	-0.0061602	3e-07	-0.0001797	-0.0024196	0.0010665	NA	-0.0000894	NA	0.0003710	NA	10	-3423.322	6866.681	31.56307	1.398533e-07
18	-2.885405	2.15e-05	-0.0063727	3e-07	0.0015667	-0.0015933	-0.0000022	NA	-0.0000495	NA	NA	-0.0001699	10	-3423.771	6867.580	32.46202	8.922136e-08
9	-2.879477	2.13e-05	-0.0065786	3e-07	-0.0000996	-0.0023764	-0.0012983	0.0e+00	NA	NA	0.0003480	NA	10	-3427.745	6875.529	40.41070	1.676618e-09
17	-2.876970	1.25e-05	-0.0065992	3e-07	0.0015377	-0.0012411	-0.0018098	1.9e-06	NA	NA	NA	-0.0002306	10	-3428.008	6876.054	40.93602	1.289324e-09
3	-2.889434	1.34e-05	-0.0061257	3e-07	0.0016199	-0.0024624	0.0007743	1.7e-06	-0.0000939	NA	NA	NA	10	-3428.372	6876.782	41.66339	8.962234e-10
12	-2.884227	2.14e-05	-0.0065599	7e-07	-0.0000916	-0.0023507	-0.0004248	NA	NA	-1e-07	0.0003405	NA	10	-3430.871	6881.780	46.66167	7.362996e-11
26	-2.893478	2.15e-05	-0.0060556	3e-07	-0.0006022	-0.0009958	0.0018498	NA	-0.0001211	NA	0.0004478	-0.0002816	11	-3430.028	6882.101	46.98312	6.269768e-11
20	-2.884117	2.15e-05	-0.0065817	7e-07	0.0015263	-0.0016648	-0.0004369	NA	NA	-1e-07	NA	-0.0001454	10	-3431.269	6882.576	47.45753	4.945782e-11
6	-2.895008	2.14e-05	-0.0061863	7e-07	0.0015803	-0.0024391	0.0018191	NA	-0.0000790	-1e-07	NA	NA	10	-3431.545	6883.128	48.01036	3.751366e-11
25	-2.877827	1.81e-05	-0.0066164	3e-07	-0.0003444	-0.0009513	-0.0015348	7.0e-07	NA	NA	0.0003863	-0.0002784	11	-3434.448	6890.942	55.82373	7.542887e-13
11	-2.891727	1.89e-05	-0.0060861	3e-07	-0.0000889	-0.0024252	0.0013324	5.0e-07	-0.0001044	NA	0.0003533	NA	11	-3434.894	6891.834	56.71575	4.828789e-13
5	-2.883424	1.60e-05	-0.0065495	7e-07	0.0015572	-0.0023903	-0.0006670	1.1e-06	NA	-1e-07	NA	NA	10	-3435.976	6891.989	56.87126	4.467550e-13
19	-2.894722	8.10e-06	-0.0058780	3e-07	0.0015805	-0.0010549	0.0020214	2.8e-06	-0.0001542	NA	NA	-0.0002814	11	-3435.056	6892.158	57.03960	4.106919e-13
28	-2.882980	2.16e-05	-0.0065992	6e-07	-0.0004288	-0.0011303	-0.0005650	NA	NA	-1e-07	0.0004000	-0.0002393	11	-3437.666	6897.377	62.25931	3.020435e-14
14	-2.901190	2.13e-05	-0.0059396	7e-07	-0.0002030	-0.0024099	0.0032037	NA	-0.0001318	-1e-07	0.0003719	NA	11	-3438.023	6898.091	62.97260	2.114363e-14
22	-2.895506	2.14e-05	-0.0061663	7e-07	0.0015517	-0.0016637	0.0019999	NA	-0.0000885	-1e-07	NA	-0.0001541	11	-3438.494	6899.034	63.91561	1.319496e-14
27	-2.898432	1.32e-05	-0.0057816	3e-07	-0.0003798	-0.0007224	0.0029159	1.7e-06	-0.0001787	NA	0.0004037	-0.0003394	12	-3441.451	6906.955	71.83674	2.513946e-16
13	-2.884235	2.15e-05	-0.0065599	7e-07	-0.0000930	-0.0023508	-0.0004230	0.0e+00	NA	-1e-07	0.0003408	NA	11	-3442.528	6907.101	71.98280	2.336900e-16
21	-2.881649	1.31e-05	-0.0065788	6e-07	0.0015150	-0.0013253	-0.0009513	1.8e-06	NA	-1e-07	NA	-0.0002092	11	-3442.806	6907.657	72.53841	1.770074e-16
7	-2.901278	1.25e-05	-0.0058625	8e-07	0.0016024	-0.0024539	0.0031392	1.8e-06	-0.0001449	-1e-07	NA	NA	11	-3443.046	6908.138	73.01974	1.391463e-16
30	-2.903249	2.15e-05	-0.0058572	7e-07	-0.0005994	-0.0010691	0.0037845	NA	-0.0001586	-1e-07	0.0004443	-0.0002653	12	-3444.762	6913.578	78.46031	9.163629e-18
29	-2.882138	1.85e-05	-0.0065971	6e-07	-0.0003202	-0.0010362	-0.0007469	7.0e-07	NA	-1e-07	0.0003770	-0.0002574	12	-3449.269	6922.592	87.47390	1.011093e-19
15	-2.903204	1.79e-05	-0.0058308	7e-07	-0.0000759	-0.0024175	0.0036222	7.0e-07	-0.0001538	-1e-07	0.0003472	NA	12	-3449.583	6923.220	88.10200	7.385839e-20
23	-2.905869	7.40e-06	-0.0056354	7e-07	0.0015654	-0.0011030	0.0042394	2.9e-06	-0.0002007	-1e-07	NA	-0.0002702	12	-3449.752	6923.557	88.43915	6.240085e-20
31	-2.909040	1.25e-05	-0.0055511	7e-07	-0.0003574	-0.0007748	0.0050239	1.9e-06	-0.0002227	-1e-07	0.0003961	-0.0003275	13	-3456.166	6938.395	103.27731	3.742176e-23

kable(dredged.ND.model, "html", caption = "ND Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

ND Dredge Table
	(Intercept)	bioclim4	log(range.size.m2)	PC1.LIG	residuals.PC1	SDi	bioclim4:SDi	log(range.size.m2):SDi	NPP:SDi	PC1.LIG:SDi	residuals.PC1:SDi	df	logLik	AICc	delta	weight
0	0.1404205	5e-07	-0.0001462	-0.0000102	-6.47e-05	-0.0000575	NA	NA	NA	NA	NA	8	14329.64	-28643.26	0.00000	9.999550e-01
8	0.1403564	4e-07	-0.0001422	-0.0000850	-5.71e-05	-0.0000618	NA	NA	NA	1.48e-05	NA	9	14320.05	-28622.06	21.19677	2.495513e-05
16	0.1404275	5e-07	-0.0001466	-0.0000108	-5.67e-05	-0.0000581	NA	NA	NA	NA	-1.6e-06	9	14319.25	-28620.46	22.79683	1.121276e-05
2	0.1400429	5e-07	-0.0001322	-0.0000107	-6.61e-05	0.0000205	NA	-2.9e-06	NA	NA	NA	9	14318.99	-28619.95	23.30593	8.692847e-06
1	0.1403637	6e-07	-0.0001457	-0.0000101	-6.47e-05	-0.0000496	0e+00	NA	NA	NA	NA	9	14314.58	-28611.13	32.12698	1.056075e-07
4	0.1406803	5e-07	-0.0001499	-0.0000090	-6.52e-05	-0.0000932	NA	NA	0	NA	NA	9	14311.61	-28605.19	38.06105	5.434109e-09
24	0.1403735	4e-07	-0.0001429	-0.0000910	-3.30e-05	-0.0000638	NA	NA	NA	1.57e-05	-4.8e-06	10	14309.74	-28599.43	43.82308	3.047333e-10
10	0.1398108	4e-07	-0.0001218	-0.0000886	-5.88e-05	0.0000503	NA	-4.2e-06	NA	1.54e-05	NA	10	14309.46	-28598.88	44.37966	2.307068e-10
18	0.1400302	5e-07	-0.0001318	-0.0000116	-5.60e-05	0.0000242	NA	-3.1e-06	NA	NA	-2.1e-06	10	14308.61	-28597.18	46.08010	9.858583e-11
9	0.1401845	8e-07	-0.0001400	-0.0000977	-5.59e-05	-0.0000401	-1e-07	NA	NA	1.74e-05	NA	10	14305.25	-28590.46	52.79550	3.432282e-12
17	0.1403708	6e-07	-0.0001459	-0.0000104	-6.13e-05	-0.0000504	0e+00	NA	NA	NA	-7.0e-07	10	14304.23	-28588.42	54.83371	1.238774e-12
3	0.1400821	5e-07	-0.0001346	-0.0000106	-6.59e-05	0.0000102	0e+00	-2.4e-06	NA	NA	NA	10	14304.00	-28587.96	55.29503	9.835981e-13
12	0.1406450	4e-07	-0.0001461	-0.0000867	-5.74e-05	-0.0001020	NA	NA	0	1.54e-05	NA	10	14302.09	-28584.14	59.11661	1.455366e-13
20	0.1406914	5e-07	-0.0001503	-0.0000097	-5.54e-05	-0.0000943	NA	NA	0	NA	-2.0e-06	10	14301.22	-28582.41	60.84729	6.125736e-14
6	0.1404437	5e-07	-0.0001414	-0.0000094	-6.61e-05	-0.0000452	NA	-1.7e-06	0	NA	NA	10	14300.95	-28581.85	61.40235	4.641190e-14
26	0.1397605	4e-07	-0.0001201	-0.0000961	-3.08e-05	0.0000624	NA	-4.8e-06	NA	1.64e-05	-5.7e-06	11	14299.18	-28576.32	66.93341	2.921291e-15
5	0.1406261	6e-07	-0.0001493	-0.0000090	-6.53e-05	-0.0000857	0e+00	NA	0	NA	NA	10	14296.55	-28573.06	70.19982	5.705361e-16
25	0.1402082	8e-07	-0.0001405	-0.0000999	-4.30e-05	-0.0000430	-1e-07	NA	NA	1.76e-05	-2.6e-06	11	14294.92	-28567.80	75.45446	4.123369e-17
11	0.1399495	7e-07	-0.0001307	-0.0000978	-5.69e-05	0.0000099	-1e-07	-2.0e-06	NA	1.73e-05	NA	11	14294.66	-28567.27	75.98128	3.168506e-17
19	0.1400559	5e-07	-0.0001333	-0.0000113	-5.77e-05	0.0000174	0e+00	-2.7e-06	NA	NA	-1.7e-06	11	14293.70	-28565.35	77.90868	1.208721e-17
28	0.1406725	4e-07	-0.0001470	-0.0000935	-3.05e-05	-0.0001055	NA	NA	0	1.64e-05	-5.4e-06	11	14291.80	-28561.55	81.70210	1.813827e-18
14	0.1402430	4e-07	-0.0001316	-0.0000890	-5.86e-05	-0.0000209	NA	-3.0e-06	0	1.58e-05	NA	11	14291.46	-28560.87	82.38244	1.290804e-18
22	0.1404322	5e-07	-0.0001410	-0.0000103	-5.50e-05	-0.0000416	NA	-1.9e-06	0	NA	-2.3e-06	11	14290.57	-28559.08	84.17121	5.277578e-19
13	0.1404724	8e-07	-0.0001438	-0.0000991	-5.62e-05	-0.0000802	-1e-07	NA	0	1.79e-05	NA	11	14287.28	-28552.52	90.73241	1.984664e-20
21	0.1406404	6e-07	-0.0001497	-0.0000094	-5.93e-05	-0.0000874	0e+00	NA	0	NA	-1.2e-06	11	14286.20	-28550.36	92.89788	6.721422e-21
7	0.1404979	6e-07	-0.0001445	-0.0000092	-6.58e-05	-0.0000590	0e+00	-1.0e-06	0	NA	NA	11	14285.96	-28549.88	93.37861	5.285327e-21
27	0.1398914	7e-07	-0.0001279	-0.0001008	-3.94e-05	0.0000252	-1e-07	-2.8e-06	NA	1.76e-05	-3.6e-06	12	14284.39	-28544.73	98.52951	4.023185e-22
30	0.1401987	4e-07	-0.0001299	-0.0000970	-2.90e-05	-0.0000096	NA	-3.5e-06	0	1.69e-05	-6.0e-06	12	14281.20	-28538.34	104.91315	1.653402e-23
29	0.1405076	8e-07	-0.0001446	-0.0001019	-4.01e-05	-0.0000847	-1e-07	NA	0	1.82e-05	-3.3e-06	12	14276.97	-28529.89	113.36490	2.416050e-25
15	0.1404182	8e-07	-0.0001418	-0.0000991	-5.65e-05	-0.0000689	-1e-07	-4.0e-07	0	1.79e-05	NA	12	14276.70	-28529.34	113.91919	1.831227e-25
23	0.1404716	5e-07	-0.0001432	-0.0000099	-5.75e-05	-0.0000518	0e+00	-1.4e-06	0	NA	-1.7e-06	12	14275.66	-28527.26	115.99255	6.494083e-26
31	0.1403601	7e-07	-0.0001390	-0.0001023	-3.86e-05	-0.0000535	-1e-07	-1.2e-06	0	1.82e-05	-3.7e-06	13	14266.43	-28506.79	136.46456	2.328506e-30

From the dredged model list we can see that the top model is one with no interactions. Run this model on the MCC and 100 trees. Noting that each model uses a unique set of values for DR/ND and a unique tree in the correlation structure.

#In both cases the top model is 1/2/3/4/5/6 no interaction terms. With no models within delta < 4: 

#Run model for DR
MCC.DR.top <- gls(MCC.DR ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(0.985, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

saveRDS(MCC.DR.top, 'data/MCC.DR.top.rds')

#Run model for ND
MCC.ND.top <- gls(MCC.ND ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corPagel(0.9996, phy = pruned.MCC.tree, fixed = TRUE), 
                data = restricted.data, 
                method = "REML")

saveRDS(MCC.ND.top, 'data/MCC.ND.top.rds')
#Run the 100 models for DR and ND using the best model:

#Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ES")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
DR.model.data <- lapply(DR.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(restricted.data$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
DR.pgls.models <- mcmapply(function(x,y) {
  gls(ES ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(0.985, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = DR.model.data, y = pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(DR.pgls.models, "data/DR.pgls.models.rds")

#Now for Node Density:
ND.model.data <- lapply(nd.list, function(x) {
  left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
ND.model.data <- lapply(ND.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
ND.pgls.models <- mcmapply(function(x,y) {
gls(ND ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(0.9996, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = ND.model.data, y = pruned.trees, 
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(ND.pgls.models, "data/ND.pgls.models.rds")

From these models we can get a distribution of the estimates for the models alongside the 95 % CIs of the MCC model, this will show the variation beween trees, and the variation associated with the top MCC model.

DR.pgls.models <- readRDS("data/DR.pgls.models.rds")
MCC.DR.top <- readRDS('data/MCC.DR.top.rds')

DR.pgls.summary <- lapply(DR.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

MCC.DR.summary <- data.frame(MCC.DR.top$coefficients, confint(MCC.DR.top)) %>% tibble::rownames_to_column()

DR.pgls.summary <- bind_rows(DR.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(MCC.DR.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `SDi` = "Sexual Dichromatism"
                    )

DR.plot <-DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

DR.plot <- DR.plot + geom_errorbarh(data = MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("DR Models")




#Now plot this for Node Density (ND):  

ND.pgls.models <- readRDS("data/ND.pgls.models.rds")
MCC.ND.top <- readRDS('data/MCC.ND.top.rds')

ND.pgls.summary <- lapply(ND.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

MCC.ND.summary <- data.frame(MCC.ND.top$coefficients, confint(MCC.ND.top)) %>% tibble::rownames_to_column()

ND.pgls.summary <- bind_rows(ND.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(MCC.ND.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

ND.plot <-ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

ND.plot <- ND.plot + geom_errorbarh(data = MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("ND Models")

PGLS Models on BAMM estimates

Using the previously generated BAMM dataframe we can now undertake model selection for speciation and extinction. The process is similar to that used on DR and ND, except given the use of a Bayesian approach in BAMM we can make use of varying levels of uncertainty between tips (species) by constructing a weighted model, where the weight is the inverse of the variance, meaning more precise estimates of speciation or extinction at a given tip (species) holds higher weight in the model.

Based on preliminary findings we found that Pagal’s lambda was = 1 and running corPagel lead to problemd of convergence. Therefore the following models are run assuming Brownian Motion with corBrownian.

MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')

#Create model dataframe for use in models
MCC.BAMM.model.data <- left_join(restricted.data %>% 
                                   dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
                                 MCC.BAMM.df %>% as.data.frame(), 
                                 by = "TipLabel")

saveRDS(MCC.BAMM.model.data, 'data/MCC.BAMM.model.data.rds')
#Prune tree
pruned.MCC.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(MCC.BAMM.model.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(MCC.BAMM.model.data) <- MCC.BAMM.model.data$TipLabel

#Run model for DR
MCC.BAMM.lambda <- gls(log(mean.lambda) ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                weights = ~ I(1/log(var.lambda)),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")

#Run model for DR
MCC.BAMM.mu <- gls(log(mean.mu) ~ SDi
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + SDi*log(range.size.m2)
                         + SDi*bioclim4
                         + SDi*residuals.PC1
                         + SDi*PC1.LIG
                         + SDi*NPP,
                weights = ~ I(1/log(var.mu)),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")


#Dredge the global MCC models

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("MCC.BAMM.model.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

dredged.MCC.lambda <- pdredge(MCC.BAMM.lambda, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(dredged.MCC.lambda, "data/dredged.MCC.lambda.rds")

dredged.MCC.mu <- pdredge(MCC.BAMM.mu, fixed = c("SDi", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(dredged.MCC.mu, "data/dredged.MCC.mu.rds")

Table S6: The dredged models for BAMM speciation and BAMM extinction both show the top model is one with no interactions , with \(\delta AICc > 4\). This is the same situation as DR/ND models (see above).

dredged.MCC.lambda <- readRDS("data/dredged.MCC.lambda.rds")
dredged.MCC.mu <- readRDS("data/dredged.MCC.mu.rds")
kable(dredged.MCC.lambda, "html", caption = "BAMM-Speciation Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

BAMM-Speciation Dredge Table
	(Intercept)	bioclim4	log(range.size.m2)	PC1.LIG	residuals.PC1	SDi	bioclim4:SDi	log(range.size.m2):SDi	NPP:SDi	PC1.LIG:SDi	residuals.PC1:SDi	df	logLik	AICc	delta	weight
0	-2.648503	-2.0e-06	-1.37e-05	-1.12e-05	3.50e-05	0.0000584	NA	NA	NA	NA	NA	8	8353.254	-16690.48	0.00000	9.999123e-01
8	-2.648548	-2.0e-06	-1.18e-05	-6.30e-05	4.19e-05	0.0000554	NA	NA	NA	1.13e-05	NA	9	8343.888	-16669.74	20.73974	3.136067e-05
16	-2.648489	-2.0e-06	-1.40e-05	-1.15e-05	4.98e-05	0.0000573	NA	NA	NA	NA	-3.0e-06	9	8343.836	-16669.64	20.84310	2.978112e-05
2	-2.647505	-2.0e-06	-5.18e-05	-8.30e-06	3.69e-05	-0.0001854	NA	9.20e-06	NA	NA	NA	9	8343.708	-16669.39	21.09789	2.621880e-05
1	-2.648473	-2.0e-06	-1.42e-05	-1.18e-05	3.51e-05	0.0000548	0e+00	NA	NA	NA	NA	9	8339.205	-16660.38	30.10536	2.901787e-07
4	-2.648288	-2.0e-06	-1.68e-05	-1.11e-05	3.85e-05	0.0000279	NA	NA	0	NA	NA	9	8336.072	-16654.11	36.36989	1.265731e-08
24	-2.648526	-2.0e-06	-1.20e-05	-7.13e-05	7.46e-05	0.0000525	NA	NA	NA	1.30e-05	-6.5e-06	10	8334.520	-16649.00	41.48242	9.821319e-10
10	-2.647644	-2.1e-06	-4.63e-05	-5.57e-05	4.30e-05	-0.0001645	NA	8.30e-06	NA	1.03e-05	NA	10	8334.336	-16648.63	41.85004	8.172261e-10
18	-2.647515	-2.0e-06	-5.12e-05	-8.50e-06	4.60e-05	-0.0001817	NA	9.00e-06	NA	NA	-1.9e-06	10	8334.293	-16648.55	41.93551	7.830375e-10
9	-2.648609	-1.9e-06	-1.06e-05	-6.59e-05	4.24e-05	0.0000621	0e+00	NA	NA	1.22e-05	NA	10	8329.896	-16639.75	50.72978	9.641194e-12
17	-2.648432	-2.1e-06	-1.50e-05	-1.25e-05	5.34e-05	0.0000505	0e+00	NA	NA	NA	-3.8e-06	10	8329.834	-16639.63	50.85452	9.058213e-12
3	-2.647431	-1.8e-06	-5.65e-05	-6.20e-06	3.71e-05	-0.0002150	0e+00	1.07e-05	NA	NA	NA	10	8329.770	-16639.50	50.98172	8.500083e-12
12	-2.648321	-2.0e-06	-1.50e-05	-6.41e-05	4.57e-05	0.0000228	NA	NA	0	1.16e-05	NA	10	8326.709	-16633.38	57.10299	3.982947e-13
20	-2.648279	-2.0e-06	-1.70e-05	-1.13e-05	5.17e-05	0.0000274	NA	NA	0	NA	-2.7e-06	10	8326.654	-16633.27	57.21322	3.769381e-13
6	-2.647016	-2.0e-06	-6.33e-05	-7.50e-06	4.20e-05	-0.0002741	NA	1.09e-05	0	NA	NA	10	8326.570	-16633.10	57.38149	3.465218e-13
26	-2.647690	-2.1e-06	-4.41e-05	-6.28e-05	6.88e-05	-0.0001512	NA	7.70e-06	NA	1.17e-05	-5.2e-06	11	8324.971	-16627.90	62.58783	2.565593e-14
5	-2.648251	-2.0e-06	-1.75e-05	-1.17e-05	3.85e-05	0.0000233	0e+00	NA	0	NA	NA	10	8322.024	-16624.01	66.47377	3.675915e-15
25	-2.648553	-2.0e-06	-1.15e-05	-7.23e-05	7.34e-05	0.0000556	0e+00	NA	NA	1.33e-05	-6.2e-06	11	8320.553	-16619.06	71.42387	3.093612e-16
11	-2.647523	-1.7e-06	-5.46e-05	-6.27e-05	4.48e-05	-0.0002202	-1e-07	1.12e-05	NA	1.28e-05	NA	11	8320.468	-16618.89	71.59252	2.843436e-16
19	-2.647437	-1.8e-06	-5.61e-05	-6.40e-06	3.94e-05	-0.0002125	0e+00	1.06e-05	NA	NA	-5.0e-07	11	8320.432	-16618.82	71.66592	2.740978e-16
28	-2.648306	-2.0e-06	-1.52e-05	-7.20e-05	7.69e-05	0.0000211	NA	NA	0	1.32e-05	-6.2e-06	11	8317.341	-16612.64	77.84794	1.245936e-17
14	-2.647149	-2.1e-06	-5.80e-05	-5.56e-05	4.83e-05	-0.0002543	NA	1.00e-05	0	1.05e-05	NA	11	8317.199	-16612.35	78.13101	1.081507e-17
22	-2.647027	-2.0e-06	-6.29e-05	-7.60e-06	4.80e-05	-0.0002708	NA	1.08e-05	0	NA	-1.2e-06	11	8317.156	-16612.27	78.21677	1.036109e-17
13	-2.648379	-1.9e-06	-1.39e-05	-6.67e-05	4.61e-05	0.0000293	0e+00	NA	0	1.24e-05	NA	11	8312.718	-16603.39	87.09376	1.224030e-19
21	-2.648217	-2.1e-06	-1.82e-05	-1.23e-05	5.56e-05	0.0000199	0e+00	NA	0	NA	-3.5e-06	11	8312.653	-16603.26	87.22383	1.146959e-19
7	-2.646914	-1.8e-06	-6.91e-05	-5.00e-06	4.24e-05	-0.0003112	0e+00	1.27e-05	0	NA	NA	11	8312.636	-16603.23	87.25674	1.128240e-19
27	-2.647565	-1.8e-06	-5.23e-05	-6.59e-05	5.94e-05	-0.0002056	-1e-07	1.05e-05	NA	1.33e-05	-3.0e-06	12	8311.150	-16598.25	92.23739	9.351196e-21
30	-2.647204	-2.1e-06	-5.56e-05	-6.18e-05	7.07e-05	-0.0002400	NA	9.50e-06	0	1.17e-05	-4.5e-06	12	8307.834	-16591.61	98.87041	3.392537e-22
29	-2.648331	-2.0e-06	-1.47e-05	-7.28e-05	7.60e-05	0.0000239	0e+00	NA	0	1.35e-05	-6.0e-06	12	8303.374	-16582.69	107.78986	3.923668e-24
15	-2.646985	-1.7e-06	-6.78e-05	-6.32e-05	5.06e-05	-0.0003213	-1e-07	1.34e-05	0	1.32e-05	NA	12	8303.340	-16582.63	107.85783	3.792547e-24
23	-2.646904	-1.8e-06	-6.96e-05	-4.80e-06	3.96e-05	-0.0003147	-1e-07	1.29e-05	0	NA	6.0e-07	12	8303.303	-16582.55	107.93102	3.656279e-24
31	-2.647020	-1.7e-06	-6.61e-05	-6.53e-05	6.00e-05	-0.0003103	-1e-07	1.29e-05	0	1.35e-05	-1.9e-06	13	8294.024	-16561.99	128.49799	1.250196e-28

kable(dredged.MCC.mu, "html", caption = "BAMM-Extinction Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

BAMM-Extinction Dredge Table
	(Intercept)	bioclim4	log(range.size.m2)	PC1.LIG	residuals.PC1	SDi	bioclim4:SDi	log(range.size.m2):SDi	NPP:SDi	PC1.LIG:SDi	residuals.PC1:SDi	df	logLik	AICc	delta	weight
0	-6.812225	4e-07	-0.0000888	-0.0001928	-1.94e-05	0.0000361	NA	NA	NA	NA	NA	8	5717.114	-11418.20	0.00000	9.998665e-01
16	-6.812223	4e-07	-0.0000891	-0.0001934	2.40e-06	0.0000343	NA	NA	NA	NA	-4.7e-06	9	5708.174	-11398.32	19.88684	4.803637e-05
8	-6.812229	4e-07	-0.0000891	-0.0001779	-2.15e-05	0.0000371	NA	NA	NA	-3.4e-06	NA	9	5708.136	-11398.24	19.96336	4.623318e-05
2	-6.811955	4e-07	-0.0000991	-0.0001921	-1.88e-05	-0.0000324	NA	2.6e-06	NA	NA	NA	9	5707.959	-11397.89	20.31734	3.873367e-05
1	-6.812299	5e-07	-0.0000876	-0.0001918	-1.94e-05	0.0000442	0e+00	NA	NA	NA	NA	9	5703.513	-11388.99	29.20883	4.542818e-07
4	-6.812053	4e-07	-0.0000912	-0.0001928	-1.74e-05	0.0000132	NA	NA	0	NA	NA	9	5700.350	-11382.67	35.53393	1.922416e-08
24	-6.812226	4e-07	-0.0000893	-0.0001827	-2.00e-06	0.0000352	NA	NA	NA	-2.4e-06	-4.1e-06	10	5699.228	-11378.42	39.78584	2.293804e-09
18	-6.811993	4e-07	-0.0000979	-0.0001928	1.60e-06	-0.0000241	NA	2.2e-06	NA	NA	-4.4e-06	10	5699.024	-11378.01	40.19459	1.869811e-09
10	-6.811924	4e-07	-0.0001008	-0.0001755	-2.11e-05	-0.0000405	NA	2.9e-06	NA	-3.8e-06	NA	10	5698.987	-11377.94	40.26855	1.801929e-09
17	-6.812267	5e-07	-0.0000884	-0.0001927	-2.00e-07	0.0000394	0e+00	NA	NA	NA	-4.1e-06	10	5694.618	-11369.20	49.00476	2.283987e-11
9	-6.812280	5e-07	-0.0000882	-0.0001801	-2.11e-05	0.0000425	0e+00	NA	NA	-2.7e-06	NA	10	5694.589	-11369.14	49.06300	2.218440e-11
3	-6.811882	6e-07	-0.0001049	-0.0001899	-1.84e-05	-0.0000705	-1e-07	4.5e-06	NA	NA	NA	10	5694.467	-11368.90	49.30697	1.963677e-11
20	-6.812053	4e-07	-0.0000915	-0.0001933	4.00e-06	0.0000117	NA	NA	0	NA	-4.6e-06	10	5691.410	-11362.78	55.42135	9.233120e-13
12	-6.812060	4e-07	-0.0000914	-0.0001784	-1.95e-05	0.0000144	NA	NA	0	-3.3e-06	NA	10	5691.372	-11362.71	55.49754	8.888011e-13
6	-6.811629	3e-07	-0.0001066	-0.0001918	-1.62e-05	-0.0000908	NA	3.7e-06	0	NA	NA	10	5691.222	-11362.41	55.79688	7.652497e-13
26	-6.811963	4e-07	-0.0000994	-0.0001801	-3.80e-06	-0.0000317	NA	2.5e-06	NA	-2.8e-06	-3.6e-06	11	5690.088	-11358.13	60.07356	9.018520e-14
5	-6.812126	5e-07	-0.0000901	-0.0001918	-1.75e-05	0.0000211	0e+00	NA	0	NA	NA	10	5686.750	-11353.46	64.74224	8.736630e-15
25	-6.812256	4e-07	-0.0000887	-0.0001837	-3.20e-06	0.0000385	0e+00	NA	NA	-2.1e-06	-3.8e-06	11	5685.710	-11349.37	68.82949	1.131900e-15
19	-6.811923	5e-07	-0.0001028	-0.0001908	-5.10e-06	-0.0000565	0e+00	3.9e-06	NA	NA	-2.9e-06	11	5685.611	-11349.18	69.02734	1.025288e-15
11	-6.811873	6e-07	-0.0001051	-0.0001791	-2.00e-05	-0.0000698	0e+00	4.5e-06	NA	-2.5e-06	NA	11	5685.544	-11349.04	69.16035	9.593193e-16
28	-6.812058	4e-07	-0.0000916	-0.0001832	-3.00e-07	0.0000127	NA	NA	0	-2.3e-06	-4.0e-06	11	5682.464	-11342.88	75.32060	4.408400e-17
22	-6.811673	4e-07	-0.0001053	-0.0001924	3.00e-06	-0.0000814	NA	3.3e-06	0	NA	-4.2e-06	11	5682.288	-11342.53	75.67391	3.694542e-17
14	-6.811598	4e-07	-0.0001083	-0.0001751	-1.85e-05	-0.0000989	NA	4.1e-06	0	-3.8e-06	NA	11	5682.250	-11342.46	75.74863	3.559064e-17
21	-6.812095	5e-07	-0.0000908	-0.0001927	1.70e-06	0.0000164	0e+00	NA	0	NA	-4.1e-06	11	5677.855	-11333.66	84.53872	4.391284e-19
13	-6.812108	5e-07	-0.0000906	-0.0001805	-1.91e-05	0.0000197	0e+00	NA	0	-2.7e-06	NA	11	5677.826	-11333.61	84.59689	4.265403e-19
7	-6.811529	6e-07	-0.0001134	-0.0001893	-1.57e-05	-0.0001358	-1e-07	6.0e-06	0	NA	NA	11	5677.734	-11333.42	84.78076	3.890758e-19
27	-6.811910	5e-07	-0.0001033	-0.0001816	-8.20e-06	-0.0000577	0e+00	3.9e-06	NA	-2.1e-06	-2.5e-06	12	5676.703	-11329.35	88.85249	5.080069e-20
30	-6.811640	4e-07	-0.0001068	-0.0001795	-2.50e-06	-0.0000895	NA	3.7e-06	0	-2.9e-06	-3.4e-06	12	5673.352	-11322.65	95.55317	1.781703e-21
29	-6.812084	4e-07	-0.0000911	-0.0001840	-1.30e-06	0.0000157	0e+00	NA	0	-2.0e-06	-3.8e-06	12	5668.947	-11313.84	104.36391	2.175735e-23
23	-6.811572	6e-07	-0.0001115	-0.0001901	-4.70e-06	-0.0001226	0e+00	5.3e-06	0	NA	-2.4e-06	12	5668.880	-11313.71	104.49749	2.035168e-23
15	-6.811524	6e-07	-0.0001135	-0.0001792	-1.72e-05	-0.0001345	0e+00	5.9e-06	0	-2.4e-06	NA	12	5668.812	-11313.57	104.63440	1.900508e-23
31	-6.811560	5e-07	-0.0001119	-0.0001812	-7.70e-06	-0.0001236	0e+00	5.4e-06	0	-2.1e-06	-2.0e-06	13	5659.972	-11293.88	124.32315	1.008124e-27

Given that we are running models weighted according to the variance of the response (\(\lambda_{BAMM}\) and \(\mu_{BAMM}\)), we can check to see whether a weighted model isw favourable to an unweighted model using an anova to compare AIC values.

#Run the top model for the MCC
#Run model for ND
MCC.Lambda.top <- gls(log(mean.lambda) ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ I(1/log(var.lambda)),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Lambda.top, 'data/MCC.Lambda.top.rds')

MCC.Mu.top <- gls(log(mean.mu) ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ I(1/log(var.mu)),
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")

saveRDS(MCC.Mu.top, 'data/MCC.Mu.top.rds')


#We can also see how the models look without the weightings: 
MCC.Lambda.top.unweighted <- gls(log(mean.lambda) ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Lambda.top.unweighted, 'data/MCC.Lambda.top.unweighted.rds')


MCC.Mu.top.unweighted <- gls(log(mean.mu) ~ SDi
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                correlation = corBrownian(phy = pruned.MCC.tree), 
                data = MCC.BAMM.model.data, 
                method = "REML")
saveRDS(MCC.Mu.top.unweighted, 'data/MCC.Mu.top.unweighted.rds')

Table S6: Models with a set of weights for the response based on the inverse of the variance of the posterior distribution in \(\lambda_{BAMM}\) and \(\mu_{BAMM}\) have much lower AIC values, indicating that the weighting scheme improves the fit of the model. In any case, we found that an unweighted model did not change the main conclusions drawn in that they did not generate any statistically significant relationships.

MCC.Lambda.top <- readRDS('data/MCC.Lambda.top.rds')
MCC.Lambda.top.unweighted <- readRDS('data/MCC.Lambda.top.unweighted.rds')
MCC.Mu.top <- readRDS('data/MCC.Mu.top.rds')
MCC.Mu.top.unweighted <- readRDS('data/MCC.Mu.top.unweighted.rds')

anova(MCC.Lambda.top, MCC.Lambda.top.unweighted) %>% pander(split.table = Inf)

	call	Model	df	AIC	BIC	logLik
MCC.Lambda.top	gls(model = log(mean.lambda) ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), weights = ~I(1/log(var.lambda)), method = “REML”)	1	8	-16691	-16637	8353
MCC.Lambda.top.unweighted	gls(model = log(mean.lambda) ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), method = “REML”)	2	8	-10674	-10620	5345

anova(MCC.Mu.top, MCC.Mu.top.unweighted) %>% pander(split.table = Inf)

	call	Model	df	AIC	BIC	logLik
MCC.Mu.top	gls(model = log(mean.mu) ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), weights = ~I(1/log(var.mu)), method = “REML”)	1	8	-11418	-11365	5717
MCC.Mu.top.unweighted	gls(model = log(mean.mu) ~ SDi + log(range.size.m2) + bioclim4 + residuals.PC1 + PC1.LIG + NPP, data = MCC.BAMM.model.data, correlation = corBrownian(phy = pruned.MCC.tree), method = “REML”)	2	8	-4024	-3971	2020

We can run the top (no interactions) model on the 100 trees, each with unique estimates of speciation and extinction from the BAMM runs.

#Read in the BAMM data for the 100 trees
BAMM.df <- readRDS('data/BAMM.df.rds')

#Take the restricted data and make it simpler with just responses and predictors.Note that we join the BAMM for the 100 trees

BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(restricted.data %>% dplyr::select(binomial, TipLabel, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
BAMM.model.data <- lapply(BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(restricted.data$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(log(mean.lambda) ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ I(1/log(var.lambda)),      
    corBrownian(phy = y), 
    data = x, 
    method = "REML")
}, x = BAMM.model.data, y = pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(BAMM.lambda.pgls.models, "data/BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(log(mean.mu) ~ SDi 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ I(1/log(var.mu)),      
    corBrownian(phy = y), 
    data = x, 
    method = "REML")
}, x = BAMM.model.data, y = pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(BAMM.mu.pgls.models, "data/BAMM.mu.pgls.models.rds")

BAMM.lambda.pgls.models <- readRDS("data/BAMM.lambda.pgls.models.rds")
MCC.Lambda.top <- readRDS('data/MCC.Lambda.top.rds')

BAMM.lambda.pgls.summary <- lapply(BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `SDi` = "Sexual Dichromatism"
                    )

MCC.lambda.summary <- data.frame(MCC.Lambda.top$coefficients, confint(MCC.Lambda.top)) %>% tibble::rownames_to_column()

BAMM.lambda.pgls.summary <- bind_rows(BAMM.lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(MCC.lambda.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

BAMM.lambda.plot <-BAMM.lambda.pgls.summary %>% filter(Parameter != "(Intercept)") %>%
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

BAMM.lambda.plot <- BAMM.lambda.plot + geom_errorbarh(data = MCC.lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 7, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("BAMM Speciation")


#Figure for extinction

BAMM.mu.pgls.models <- readRDS("data/BAMM.mu.pgls.models.rds")
MCC.Mu.top <- readRDS('data/MCC.Mu.top.rds')

BAMM.mu.pgls.summary <- lapply(BAMM.mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

MCC.mu.summary <- data.frame(MCC.Mu.top$coefficients, confint(MCC.Mu.top)) %>% tibble::rownames_to_column()

BAMM.mu.pgls.summary <- bind_rows(BAMM.mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(MCC.mu.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

BAMM.mu.plot <-BAMM.mu.pgls.summary %>% filter(Parameter != "(Intercept)") %>%
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

BAMM.mu.plot <- BAMM.mu.plot + geom_errorbarh(data = MCC.mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = MCC.mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,0.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 7, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,0.5,.5),"cm"))+
  ggtitle("BAMM Extinction")

PGLS Model Summary

grid.arrange(
  symmetrise_scale(DR.plot, "x"),
  symmetrise_scale(ND.plot, "x"),
  symmetrise_scale(BAMM.lambda.plot, "x"),
  symmetrise_scale(BAMM.mu.plot, "x"),
  nrow = 1
)

Figure S8: This figure is the same basic figure as seen in the manuscrip (FIgure 1), it provides model estimates for four response variables across 100 random trees alongside the MCC tree.

From the estimates of 100 trees we can calculate HPD intervals. Note that these do no account for the CI of each estimate. Thus it is a more a representation of tree uncertainty than model uncertainty.

Table S7: Highest posterior density intervals are calculated for the above figure (ie; models using RGB values of sexual dichromatism). The HPD range is determined using the hdi function with a 95 % credible interval. These intervals do not take into account the variance associated with each interval and thus are not an estimate of model precision. Intervals not overlapping zero suggest that 95 % of trees from the posterior generate a model estimate for the given parameter that in the same direction (+ or -).

hpd.DR.top <- list()
for(x in unique(DR.pgls.summary$Parameter)) {
hpd.DR.top[[x]] = hdi(DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select("Estimate"))
}
hpd.DR.top <- bind_rows(hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)")


saveRDS(hpd.DR.top, 'data/hpd.DR.top.rds')

#For ND
hpd.ND.top <- list()
for(x in unique(ND.pgls.summary$Parameter)) {
hpd.ND.top[[x]] = hdi(ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.ND.top <- bind_rows(hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)") 

saveRDS(hpd.ND.top, 'data/hpd.ND.top.rds')

hpd.Lambda.top <- list()
for(x in unique(BAMM.lambda.pgls.summary$Parameter)) {
hpd.Lambda.top[[x]] = hdi(BAMM.lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.Lambda.top <- bind_rows(hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)") 


saveRDS(hpd.Lambda.top, 'data/hpd.Lambda.top.rds')

hpd.Mu.top <- list()
for(x in unique(BAMM.mu.pgls.summary$Parameter)) {
hpd.Mu.top[[x]] = hdi(BAMM.mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
hpd.Mu.top <- bind_rows(hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper")) %>% dplyr::select(-"(Intercept)")


saveRDS(hpd.Mu.top, 'data/hpd.Mu.top.rds')

hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "DR HPD Interval")

DR HPD Interval
	SDi	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.00193	-0.00844	-2.16e-05	-0.00323	-0.00389	-1.4e-06
Upper	0.00156	-0.00182	5.6e-05	0.00479	0.00523	1.53e-06

hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "ND HPD Interval")

ND HPD Interval
	SDi	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-7.68e-05	-0.00019	-9.72e-07	-9.09e-05	-0.00015	-4.27e-08
Upper	5.97e-05	9.11e-06	1.62e-06	0.000175	0.000117	4.28e-08

hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "BAMM Speciation HPD Interval")

BAMM Speciation HPD Interval
	SDi	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.00637	-0.00453	-0.000144	-0.00606	-0.0109	-3.18e-06
Upper	0.0037	0.00856	5.74e-05	0.0136	0.00699	2.99e-06

hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "BAMM Extinction HPD Interval")

BAMM Extinction HPD Interval
	SDi	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.00539	-0.00747	-0.000134	-0.0137	-0.0118	-3.05e-06
Upper	0.00531	0.0124	0.000122	0.0155	0.0114	3.18e-06

Subsetted analysis with Armenta data

Using the dataset from Armenta et al. (2008) we can conduct a subsetted analysis where dichromatism values are measured from spectrophotmetry and synthesised into a colour discriminability measure, which is a difference measure between the sexes. To make this dataset comprable with ours we use the absolute difference between the sexes. Therefore making the scale from monochromatism to dichromatism rather than female colouration to male colouration.

#Read in Armenta data
Armenta.data <- read.csv('data/Armenta_2008.csv', stringsAsFactors = F)
MCC.BAMM.df <- readRDS('data/MCC.BAMM.df.rds')
#A couple of rows have "-", remove from dataset
Armenta.data <- Armenta.data %>% dplyr::select(binomial, Colour.discriminability) %>% mutate(Colour.discriminability = replace(Colour.discriminability, Colour.discriminability == "-", "NA")) %>% filter(Colour.discriminability != "NA")
Armenta.data$Colour.discriminability <- Armenta.data$Colour.discriminability %>% as.numeric()

MCC.Armenta.model.data <- inner_join(restricted.data %>% dplyr::select(TipLabel, binomial, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP, MCC.DR, MCC.ND), Armenta.data, by = "binomial")

MCC.Armenta.model.data <- inner_join(MCC.Armenta.model.data, MCC.BAMM.df, by = "TipLabel") %>% filter(Colour.discriminability != "NA")
MCC.Armenta.model.data$abs.Colour.discriminability <- abs(MCC.Armenta.model.data$Colour.discriminability)

#Plot the correlation
p1 <- MCC.Armenta.model.data %>% ggplot(aes(x = SDi, y = abs(Colour.discriminability)))+
  geom_point()+
  geom_smooth(method = 'loess')+
  theme_minimal()

grid.newpage()
ggExtra::ggMarginal(
  p = p1,
  type = 'density',
  margins = 'both',
  size = 5,
  colour = 'black',
  fill = 'gray'
)

Figure S9: There is a correlation between the RGB measures and the Spectophotometry measures, The RGB measures seem to be more noisy around the lower values, where s

Run PGLS model for the data

#Prune tree
pruned.MCC.Armenta.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(MCC.Armenta.model.data$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(MCC.Armenta.model.data) <- MCC.Armenta.model.data$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.Armenta.global <- gls(MCC.DR ~ abs.Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + abs.Colour.discriminability*log(range.size.m2)
                         + abs.Colour.discriminability*bioclim4
                         + abs.Colour.discriminability*residuals.PC1
                         + abs.Colour.discriminability*PC1.LIG
                         + abs.Colour.discriminability*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = FALSE), 
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.ND.Armenta.global <- gls(MCC.ND ~ abs.Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + abs.Colour.discriminability*log(range.size.m2)
                         + abs.Colour.discriminability*bioclim4
                         + abs.Colour.discriminability*residuals.PC1
                         + abs.Colour.discriminability*PC1.LIG
                         + abs.Colour.discriminability*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), #lambda = 1
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Lambda.Armenta.global <- gls(log(mean.lambda) ~ abs.Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + abs.Colour.discriminability*log(range.size.m2)
                         + abs.Colour.discriminability*bioclim4
                         + abs.Colour.discriminability*residuals.PC1
                         + abs.Colour.discriminability*PC1.LIG
                         + abs.Colour.discriminability*NPP,
                weights = ~ I(1/log(var.lambda)),
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), 
                data = MCC.Armenta.model.data, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Extinction.Armenta.global <- gls(log(mean.mu) ~ abs.Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + abs.Colour.discriminability*log(range.size.m2)
                         + abs.Colour.discriminability*bioclim4
                         + abs.Colour.discriminability*residuals.PC1
                         + abs.Colour.discriminability*PC1.LIG
                         + abs.Colour.discriminability*NPP,
                weights = ~ I(1/log(var.mu)),
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), 
                data = MCC.Armenta.model.data, 
                method = "REML")

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("MCC.Armenta.model.data"), envir=environment())
clusterExport(myclust, c("pruned.MCC.Armenta.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

Armenta.dredged.ND.model <- pdredge(MCC.ND.Armenta.global, fixed = c("abs.Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.ND.model, 'data/Armenta.dredged.ND.model.rds')

Armenta.dredged.DR.model <- pdredge(MCC.DR.Armenta.global, fixed = c("abs.Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.DR.model, 'data/Armenta.dredged.DR.model.rds')

Armenta.dredged.spec.model <- pdredge(MCC.Lambda.Armenta.global, fixed = c("abs.Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.spec.model, 'data/Armenta.dredged.spec.model.rds')

Armenta.dredged.extinct.model <- pdredge(MCC.Extinction.Armenta.global, fixed = c("abs.Colour.discriminability", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Armenta.dredged.extinct.model, 'data/Armenta.dredged.extinct.model.rds')

Table S8: The dredged models for Models using a subset dataset where sexual dichromatism was measured using spectrophotometry, all four show the top model is one with no interactions , with \(\delta AICc > 4\). This is the same situation as DR/ND models (see above).

Armenta.dredged.ND.model <- readRDS("data/Armenta.dredged.ND.model.rds")
Armenta.dredged.DR.model <- readRDS("data/Armenta.dredged.DR.model.rds")
Armenta.dredged.spec.model <- readRDS("data/Armenta.dredged.spec.model.rds")
Armenta.dredged.extinct.model <- readRDS("data/Armenta.dredged.extinct.model.rds")


kable(Armenta.dredged.ND.model, "html", caption = "ND Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

ND Spectrophotometry Dichromatism Dredge Table
	(Intercept)	abs.Colour.discriminability	bioclim4	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	abs.Colour.discriminability:bioclim4	abs.Colour.discriminability:log(range.size.m2)	abs.Colour.discriminability:NPP	abs.Colour.discriminability:PC1.LIG	abs.Colour.discriminability:residuals.PC1	df	logLik	AICc	delta	weight
0	0.2077812	-0.0002188	7.0e-06	-0.0000307	5e-07	-0.0008917	0.0010513	NA	NA	NA	NA	NA	8	1033.5047	-2050.758	0.00000	9.974822e-01
2	0.2021525	0.0082305	6.6e-06	0.0001692	5e-07	-0.0008672	0.0010498	NA	-0.0002958	NA	NA	NA	9	1027.5845	-2036.854	13.90384	9.543858e-04
8	0.2080321	-0.0000682	6.7e-06	-0.0000426	5e-07	-0.0005900	0.0010626	NA	NA	NA	-0.0004021	NA	9	1027.4340	-2036.553	14.20484	8.210378e-04
16	0.2072009	-0.0003166	7.0e-06	-0.0000055	5e-07	-0.0009221	0.0007620	NA	NA	NA	NA	0.0003857	9	1027.3188	-2036.322	14.43529	7.316795e-04
1	0.2085571	-0.0024081	4.3e-06	-0.0000037	5e-07	-0.0009627	0.0010369	4.80e-06	NA	NA	NA	NA	9	1022.7301	-2027.145	23.61265	7.438452e-06
24	0.2071343	-0.0000863	6.6e-06	-0.0000048	5e-07	-0.0002820	0.0004750	NA	NA	NA	-0.0008968	0.0008019	10	1021.7263	-2023.067	27.69092	9.680505e-07
10	0.2064178	0.0023155	6.7e-06	0.0000149	5e-07	-0.0006057	0.0010613	NA	-0.0000838	NA	-0.0003719	NA	10	1021.5903	-2022.795	27.96290	8.449648e-07
18	0.2011316	0.0087725	6.7e-06	0.0002101	5e-07	-0.0008964	0.0007539	NA	-0.0003183	NA	NA	0.0003943	10	1021.4064	-2022.427	28.33081	7.029863e-07
4	0.2090202	-0.0039708	6.6e-06	0.0000038	2e-07	-0.0008146	0.0010507	NA	NA	4e-07	NA	NA	9	1020.3173	-2022.319	28.43836	6.661846e-07
9	0.2119186	-0.0072833	-3.8e-06	0.0000024	4e-07	0.0004728	0.0010599	1.73e-05	NA	NA	-0.0021588	NA	10	1018.4257	-2016.465	34.29215	3.568279e-08
3	0.1906298	0.0238655	1.6e-06	0.0006668	4e-07	-0.0009267	0.0010228	7.80e-06	-0.0009674	NA	NA	NA	10	1017.2485	-2014.111	36.64656	1.099525e-08
17	0.2082477	-0.0021612	4.6e-06	0.0000019	5e-07	-0.0009645	0.0009343	4.20e-06	NA	NA	NA	0.0001392	10	1016.5624	-2012.739	38.01879	5.536429e-09
26	0.2115099	-0.0066228	6.8e-06	-0.0001604	5e-07	-0.0002217	0.0004455	NA	0.0002299	NA	-0.0010072	0.0008469	11	1015.9362	-2009.408	41.34925	1.047198e-09
6	0.2116837	-0.0081513	6.7e-06	-0.0000832	2e-07	-0.0008182	0.0010514	NA	0.0001337	5e-07	NA	NA	10	1014.5280	-2008.670	42.08750	7.239709e-10
12	0.2090246	-0.0034064	6.5e-06	-0.0000081	3e-07	-0.0006328	0.0010580	NA	NA	4e-07	-0.0002551	NA	10	1014.2134	-2008.041	42.71676	5.285413e-10
20	0.2084568	-0.0043131	6.6e-06	0.0000338	2e-07	-0.0008433	0.0007312	NA	NA	4e-07	NA	0.0004259	10	1014.1659	-2007.946	42.81174	5.040264e-10
11	0.1976270	0.0136600	-5.6e-06	0.0005318	4e-07	0.0004453	0.0010479	1.92e-05	-0.0007642	NA	-0.0020748	NA	11	1012.8242	-2003.184	47.57326	4.661264e-11
25	0.2110580	-0.0068143	-3.2e-06	0.0000248	4e-07	0.0006083	0.0006664	1.61e-05	NA	NA	-0.0023731	0.0005373	11	1012.4755	-2002.487	48.27063	3.289058e-11
19	0.1902052	0.0245924	1.3e-06	0.0006857	4e-07	-0.0009247	0.0010743	8.20e-06	-0.0009988	NA	NA	-0.0000705	11	1011.1192	-1999.774	50.98317	8.473279e-12
5	0.2112711	-0.0105178	1.7e-06	0.0000755	1e-07	-0.0008747	0.0010260	8.10e-06	NA	7e-07	NA	NA	10	1010.0432	-1999.700	51.05717	8.165496e-12
14	0.2152129	-0.0129387	6.7e-06	-0.0002139	2e-07	-0.0005837	0.0010617	NA	0.0003106	4e-07	-0.0003358	NA	11	1008.5207	-1994.577	56.18019	6.302792e-13
28	0.2080277	-0.0029321	6.4e-06	0.0000227	3e-07	-0.0003346	0.0005018	NA	NA	3e-07	-0.0007456	0.0007600	11	1008.4601	-1994.456	56.30141	5.932118e-13
22	0.2111255	-0.0085019	6.7e-06	-0.0000534	2e-07	-0.0008469	0.0007319	NA	0.0001340	5e-07	NA	0.0004259	11	1008.3773	-1994.291	56.46701	5.460721e-13
13	0.2150226	-0.0163349	-6.9e-06	0.0000891	0e+00	0.0006307	0.0010490	2.14e-05	NA	8e-07	-0.0022521	NA	11	1005.8636	-1989.263	61.49446	4.421320e-14
27	0.2006302	0.0086723	-4.7e-06	0.0004135	4e-07	0.0005527	0.0007595	1.78e-05	-0.0005695	NA	-0.0022549	0.0003979	12	1006.8353	-1989.121	61.63628	4.118677e-14
7	0.2018300	0.0039116	8.0e-07	0.0004027	1e-07	-0.0008692	0.0010204	9.10e-06	-0.0004885	6e-07	NA	NA	11	1004.3954	-1986.327	64.43070	1.018491e-14
21	0.2113581	-0.0106230	1.6e-06	0.0000747	1e-07	-0.0008738	0.0010493	8.20e-06	NA	7e-07	NA	-0.0000318	11	1003.8780	-1985.292	65.46550	6.070883e-15
30	0.2208911	-0.0228961	6.8e-06	-0.0004056	2e-07	-0.0001895	0.0004313	NA	0.0006526	5e-07	-0.0009842	0.0008670	12	1002.8892	-1981.229	69.52846	7.961453e-16
15	0.2117492	-0.0113312	-7.2e-06	0.0002010	0e+00	0.0006174	0.0010469	2.16e-05	-0.0001674	8e-07	-0.0022294	NA	12	1000.1626	-1975.776	74.98178	5.209781e-17
29	0.2142549	-0.0154870	-6.3e-06	0.0000992	0e+00	0.0007138	0.0007815	2.04e-05	NA	8e-07	-0.0023927	0.0003661	12	999.8185	-1975.088	75.66996	3.693025e-17
23	0.2012141	0.0049585	3.0e-07	0.0004324	1e-07	-0.0008651	0.0011085	9.80e-06	-0.0005374	6e-07	NA	-0.0001206	12	998.2747	-1972.000	78.75754	7.887190e-18
31	0.2143554	-0.0156416	-6.3e-06	0.0000957	0e+00	0.0007145	0.0007808	2.03e-05	0.0000053	8e-07	-0.0023938	0.0003672	13	994.1598	-1961.678	89.08005	4.522906e-20

kable(Armenta.dredged.DR.model, "html", caption = "DR Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

DR Spectrophotometry Dichromatism Dredge Table
	(Intercept)	abs.Colour.discriminability	bioclim4	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	abs.Colour.discriminability:bioclim4	abs.Colour.discriminability:log(range.size.m2)	abs.Colour.discriminability:NPP	abs.Colour.discriminability:PC1.LIG	abs.Colour.discriminability:residuals.PC1	df	logLik	AICc	delta	weight
0	-2.686495	0.0070877	0.0001787	0.0092136	7.8e-06	-0.0316161	0.0155979	NA	NA	NA	NA	NA	9	-512.1586	1042.632	0.000000	9.381493e-01
2	-3.007955	0.4954768	0.0001616	0.0207305	7.4e-06	-0.0303962	0.0155527	NA	-0.0171356	NA	NA	NA	10	-514.8054	1049.997	7.364384	2.361111e-02
8	-2.671937	0.0095478	0.0001664	0.0088630	7.5e-06	-0.0204339	0.0165739	NA	NA	NA	-0.0141963	NA	10	-514.9950	1050.376	7.743674	1.953235e-02
16	-2.693478	0.0026019	0.0001826	0.0096006	7.7e-06	-0.0329197	0.0067574	NA	NA	NA	NA	0.0119257	10	-515.1631	1050.712	8.079823	1.651052e-02
24	-2.671555	0.0021899	0.0001618	0.0093363	6.9e-06	-0.0108082	-0.0024601	NA	NA	NA	-0.0302019	0.0271157	11	-516.7325	1055.929	13.296582	1.216062e-03
18	-3.030254	0.5136141	0.0001649	0.0216731	7.3e-06	-0.0317065	0.0062567	NA	-0.0179362	NA	NA	0.0125392	11	-517.7659	1057.996	15.363318	4.326830e-04
10	-2.888922	0.3319758	0.0001586	0.0165913	7.3e-06	-0.0229029	0.0162654	NA	-0.0113381	NA	-0.0100322	NA	11	-517.8235	1058.111	15.478518	4.084646e-04
1	-2.684460	-0.0005100	0.0001701	0.0092809	7.8e-06	-0.0319468	0.0154800	0.0000173	NA	NA	NA	NA	10	-520.3848	1061.156	18.523187	8.912774e-05
26	-2.708387	0.0570878	0.0001604	0.0106474	6.9e-06	-0.0113655	-0.0022010	NA	-0.0019260	NA	-0.0292503	0.0266915	12	-519.6806	1063.911	21.278224	2.247835e-05
4	-2.664734	-0.0766697	0.0001737	0.0101956	2.5e-06	-0.0302954	0.0151253	NA	NA	9.40e-06	NA	NA	10	-522.3115	1065.009	22.376701	1.297876e-05
9	-2.620545	-0.0880714	0.0000330	0.0091232	7.6e-06	-0.0063902	0.0165039	0.0002312	NA	NA	-0.0377099	NA	11	-521.8047	1066.073	23.440942	7.623179e-06
3	-3.198180	0.7554332	0.0000956	0.0284381	7.4e-06	-0.0317911	0.0147218	0.0001126	-0.0279927	NA	NA	NA	11	-522.4668	1067.398	24.765264	3.931546e-06
17	-2.698378	0.0162276	0.0001990	0.0095441	7.6e-06	-0.0324836	0.0056685	-0.0000325	NA	NA	NA	0.0137191	11	-523.2506	1068.965	26.332758	1.795504e-06
25	-2.624345	-0.0878830	0.0000385	0.0095694	7.1e-06	0.0017317	-0.0015856	0.0002142	NA	NA	-0.0511968	0.0257625	12	-523.6822	1071.914	29.281408	4.110503e-07
11	-3.103273	0.6223668	-0.0000320	0.0270494	7.3e-06	-0.0071098	0.0157656	0.0003129	-0.0261909	NA	-0.0364040	NA	12	-523.9878	1072.525	29.892508	3.028283e-07
6	-2.872767	0.2560291	0.0001646	0.0170736	4.0e-06	-0.0299649	0.0152592	NA	-0.0107060	6.30e-06	NA	NA	11	-525.1190	1072.702	30.069689	2.771545e-07
12	-2.656835	-0.0592550	0.0001645	0.0097310	3.2e-06	-0.0213688	0.0160145	NA	NA	7.70e-06	-0.0116416	NA	11	-525.3025	1073.069	30.436578	2.307029e-07
20	-2.672292	-0.0770310	0.0001777	0.0105126	2.6e-06	-0.0315844	0.0067948	NA	NA	9.00e-06	NA	0.0112662	11	-525.3603	1073.185	30.552288	2.177344e-07
19	-3.136655	0.6624291	0.0001249	0.0259779	7.3e-06	-0.0321603	0.0083676	0.0000666	-0.0241379	NA	NA	0.0090409	12	-525.5520	1075.653	33.021033	6.336447e-08
28	-2.663124	-0.0357086	0.0001610	0.0097929	4.6e-06	-0.0118446	-0.0017365	NA	NA	4.30e-06	-0.0279130	0.0256460	12	-527.2371	1079.024	36.391183	1.174970e-08
27	-2.923316	0.3530486	-0.0000027	0.0206285	6.9e-06	0.0000782	0.0006450	0.0002671	-0.0162510	NA	-0.0483668	0.0219478	13	-526.2886	1079.219	36.586890	1.065442e-08
22	-2.914834	0.3101861	0.0001673	0.0185352	4.4e-06	-0.0312776	0.0064477	NA	-0.0124597	5.40e-06	NA	0.0119474	12	-528.1204	1080.790	38.157781	4.857527e-09
14	-2.743795	0.0765554	0.0001616	0.0126339	3.8e-06	-0.0221897	0.0159845	NA	-0.0044300	6.60e-06	-0.0103857	NA	12	-528.1225	1080.794	38.161874	4.847595e-09
5	-2.646902	-0.1435694	0.0001311	0.0108597	8.0e-07	-0.0314027	0.0143755	0.0000830	NA	1.28e-05	NA	NA	11	-530.1941	1082.852	40.219766	1.732450e-09
13	-2.555957	-0.3017598	-0.0000487	0.0112730	-2.1e-06	-0.0003634	0.0152499	0.0003642	NA	1.76e-05	-0.0453855	NA	12	-531.0085	1086.566	43.933901	2.704871e-10
30	-2.603416	-0.1287498	0.0001626	0.0078108	4.2e-06	-0.0111079	-0.0019637	NA	0.0030497	5.00e-06	-0.0290449	0.0260425	13	-530.0513	1086.745	44.112144	2.474239e-10
7	-3.045354	0.4701151	0.0000871	0.0247962	2.7e-06	-0.0314447	0.0141569	0.0001342	-0.0210759	8.90e-06	NA	NA	12	-532.6233	1089.796	47.163630	5.380466e-11
21	-2.662261	-0.1088665	0.0001564	0.0107721	1.8e-06	-0.0318399	0.0081816	0.0000396	NA	1.07e-05	NA	0.0089216	12	-533.2753	1091.100	48.467514	2.803399e-11
29	-2.573270	-0.2544878	-0.0000263	0.0111512	-4.0e-07	0.0047402	0.0011884	0.0003213	NA	1.37e-05	-0.0543662	0.0204096	13	-533.4178	1093.478	50.845155	8.538604e-12
15	-2.831379	0.1234432	-0.0000707	0.0207451	-7.0e-07	-0.0017671	0.0150646	0.0003865	-0.0143568	1.47e-05	-0.0433697	NA	13	-533.6497	1093.941	51.309030	6.771077e-12
23	-3.021846	0.4420754	0.0001096	0.0235154	3.3e-06	-0.0317586	0.0096218	0.0000977	-0.0192309	7.60e-06	NA	0.0065712	13	-535.7653	1098.173	55.540280	8.163068e-13
31	-2.716186	-0.0349812	-0.0000391	0.0161109	3.0e-07	0.0036812	0.0019748	0.0003353	-0.0075020	1.24e-05	-0.0527335	0.0191444	14	-536.1666	1101.075	58.442891	1.912315e-13

kable(Armenta.dredged.spec.model, "html", caption = "BAMM Speciation Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

BAMM Speciation Spectrophotometry Dichromatism Dredge Table
	(Intercept)	abs.Colour.discriminability	bioclim4	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	abs.Colour.discriminability:bioclim4	abs.Colour.discriminability:log(range.size.m2)	abs.Colour.discriminability:NPP	abs.Colour.discriminability:PC1.LIG	abs.Colour.discriminability:residuals.PC1	df	logLik	AICc	delta	weight
0	-3.102033	0.0005342	-3.25e-05	0.0016996	-1.4e-06	-0.0019005	0.0021590	NA	NA	NA	NA	NA	8	254.4471	-492.6424	0.000000	9.865856e-01
2	-3.011015	-0.1423896	-2.74e-05	-0.0016071	-1.3e-06	-0.0023624	0.0022014	NA	0.0049897	NA	NA	NA	9	250.6065	-482.8977	9.744703	7.552653e-03
8	-3.102034	0.0005328	-3.25e-05	0.0016996	-1.4e-06	-0.0019034	0.0021588	NA	NA	NA	0.0000038	NA	9	249.6525	-480.9898	11.652618	2.909383e-03
16	-3.100824	0.0011510	-3.30e-05	0.0016302	-1.4e-06	-0.0017539	0.0032438	NA	NA	NA	NA	-0.0014356	9	249.6368	-480.9583	11.684104	2.863940e-03
10	-2.984592	-0.1824367	-2.76e-05	-0.0025773	-1.3e-06	-0.0006878	0.0023634	NA	0.0064189	NA	-0.0024181	NA	10	246.1462	-471.9064	20.735999	3.100059e-05
18	-3.007807	-0.1446367	-2.77e-05	-0.0017519	-1.3e-06	-0.0022091	0.0034060	NA	0.0050920	NA	NA	-0.0015931	10	245.8256	-471.2652	21.377207	2.249745e-05
1	-3.102052	0.0006060	-3.25e-05	0.0016988	-1.4e-06	-0.0018980	0.0021598	-2.00e-07	NA	NA	NA	NA	9	244.7581	-471.2009	21.441484	2.178591e-05
24	-3.100506	0.0009393	-3.22e-05	0.0016134	-1.3e-06	-0.0026428	0.0036166	NA	NA	NA	0.0012712	-0.0020332	10	245.0686	-469.7511	22.891224	1.055282e-05
4	-3.101388	-0.0013224	-3.27e-05	0.0017160	-1.5e-06	-0.0018653	0.0021526	NA	NA	2.0e-07	NA	NA	9	242.3609	-466.4066	26.235779	1.982019e-06
3	-2.975150	-0.1922942	-1.26e-05	-0.0031422	-1.3e-06	-0.0021749	0.0023320	-2.33e-05	0.0071143	NA	NA	NA	10	241.4066	-462.4272	30.215135	2.710202e-07
26	-2.989499	-0.1737415	-2.78e-05	-0.0024101	-1.3e-06	-0.0010229	0.0029027	NA	0.0061192	NA	-0.0018279	-0.0007655	11	241.4495	-460.4351	32.207318	1.000932e-07
9	-3.102090	0.0007097	-3.23e-05	0.0016978	-1.4e-06	-0.0019276	0.0021582	-4.00e-07	NA	NA	0.0000452	NA	10	240.3903	-460.3947	32.247714	9.809181e-08
17	-3.099659	-0.0019395	-3.68e-05	0.0016502	-1.4e-06	-0.0018305	0.0034866	6.90e-06	NA	NA	NA	-0.0018010	10	240.0956	-459.8052	32.837151	7.305301e-08
6	-2.965005	-0.2236032	-2.65e-05	-0.0028552	-2.4e-06	-0.0021865	0.0021511	NA	0.0071363	2.2e-06	NA	NA	10	239.0351	-457.6843	34.958039	2.529840e-08
12	-3.101362	-0.0014584	-3.26e-05	0.0017177	-1.5e-06	-0.0019218	0.0021473	NA	NA	2.0e-07	0.0000783	NA	10	237.5941	-454.8022	37.840122	5.987653e-09
20	-3.100049	-0.0010492	-3.31e-05	0.0016492	-1.5e-06	-0.0017109	0.0032446	NA	NA	2.0e-07	NA	-0.0014467	10	237.5537	-454.7213	37.921034	5.750252e-09
11	-2.973061	-0.1962325	-1.53e-05	-0.0031810	-1.3e-06	-0.0016808	0.0023603	-1.93e-05	0.0071952	NA	-0.0007606	NA	11	237.0550	-451.6461	40.996266	1.235687e-09
19	-2.978672	-0.1865916	-1.47e-05	-0.0029942	-1.3e-06	-0.0021440	0.0027509	-2.03e-05	0.0068739	NA	NA	-0.0005765	11	236.6406	-450.8172	41.825217	8.164026e-10
25	-3.100294	0.0003062	-3.31e-05	0.0016195	-1.3e-06	-0.0025599	0.0036266	1.50e-06	NA	NA	0.0011294	-0.0020442	11	235.8111	-449.1582	43.484205	3.561719e-10
14	-2.936793	-0.2665863	-2.68e-05	-0.0038823	-2.5e-06	-0.0004524	0.0023176	NA	0.0086555	2.3e-06	-0.0024991	NA	11	234.5894	-446.7148	45.927584	1.049752e-10
22	-2.959170	-0.2301133	-2.69e-05	-0.0030776	-2.5e-06	-0.0020078	0.0034852	NA	0.0073560	2.3e-06	NA	-0.0017674	11	234.2914	-446.1189	46.523451	7.792839e-11
5	-3.101147	-0.0021301	-3.32e-05	0.0017245	-1.5e-06	-0.0018753	0.0021466	1.00e-06	NA	2.0e-07	NA	NA	10	232.7712	-445.1564	47.486009	4.815916e-11
28	-3.098868	-0.0035688	-3.24e-05	0.0016485	-1.6e-06	-0.0027373	0.0036944	NA	NA	5.0e-07	0.0015312	-0.0021780	11	233.0483	-443.6326	49.009797	2.247980e-11
27	-2.976536	-0.1904124	-1.60e-05	-0.0030519	-1.3e-06	-0.0018275	0.0026779	-1.83e-05	0.0069793	NA	-0.0004974	-0.0004506	12	232.3627	-440.1761	52.466258	3.992378e-12
7	-2.942202	-0.2527996	-1.45e-05	-0.0039375	-2.2e-06	-0.0020578	0.0022663	-1.92e-05	0.0085725	1.9e-06	NA	NA	11	229.7460	-437.0281	55.614256	8.272802e-13
30	-2.942035	-0.2572896	-2.69e-05	-0.0036970	-2.5e-06	-0.0008737	0.0030026	NA	0.0083212	2.3e-06	-0.0017500	-0.0009736	12	229.9093	-435.2692	57.373161	3.433295e-13
13	-3.101109	-0.0022357	-3.34e-05	0.0017254	-1.5e-06	-0.0018530	0.0021476	1.20e-06	NA	3.0e-07	-0.0000338	NA	11	228.4163	-434.3686	58.273768	2.188501e-13
21	-3.096804	-0.0099118	-3.95e-05	0.0017149	-1.7e-06	-0.0017578	0.0036361	1.11e-05	NA	7.0e-07	NA	-0.0020538	11	228.1812	-433.8984	58.743975	1.729986e-13
15	-2.934658	-0.2668675	-2.02e-05	-0.0040951	-2.4e-06	-0.0010031	0.0023198	-1.03e-05	0.0088820	2.1e-06	-0.0016061	NA	12	225.4547	-426.3600	66.282335	3.991248e-15
23	-2.945869	-0.2475112	-1.85e-05	-0.0037298	-2.3e-06	-0.0019886	0.0030527	-1.31e-05	0.0082503	2.1e-06	NA	-0.0010905	12	225.0366	-425.5238	67.118540	2.627417e-15
29	-3.097523	-0.0074677	-3.60e-05	0.0016836	-1.7e-06	-0.0024325	0.0037579	5.90e-06	NA	7.0e-07	0.0010393	-0.0022656	12	223.8961	-423.2430	69.399375	8.399473e-16
31	-2.939504	-0.2589943	-2.17e-05	-0.0038954	-2.4e-06	-0.0012438	0.0028952	-8.20e-06	0.0085543	2.2e-06	-0.0011604	-0.0008185	13	220.7883	-414.9347	77.707705	1.318624e-17

kable(Armenta.dredged.extinct.model, "html", caption = "BAMM Extinction Spectrophotometry Dichromatism Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

BAMM Extinction Spectrophotometry Dichromatism Dredge Table
	(Intercept)	abs.Colour.discriminability	bioclim4	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	abs.Colour.discriminability:bioclim4	abs.Colour.discriminability:log(range.size.m2)	abs.Colour.discriminability:NPP	abs.Colour.discriminability:PC1.LIG	abs.Colour.discriminability:residuals.PC1	df	logLik	AICc	delta	weight
0	-6.792943	-0.0103858	3.41e-05	-0.0002965	-8.0e-07	-0.0106156	0.0000352	NA	NA	NA	NA	NA	8	-80.41697	177.0857	0.000000	9.755486e-01
2	-6.629991	-0.2696903	4.33e-05	-0.0062547	-6.0e-07	-0.0114938	0.0000244	NA	0.0090568	NA	NA	NA	9	-83.65433	185.6239	8.538207	1.365214e-02
16	-6.791609	-0.0089665	3.27e-05	-0.0004075	-8.0e-07	-0.0102634	0.0025209	NA	NA	NA	NA	-0.0033232	9	-84.56739	187.4500	10.364333	5.478512e-03
8	-6.792929	-0.0103492	3.41e-05	-0.0002974	-8.0e-07	-0.0105457	0.0000407	NA	NA	NA	-0.0000938	NA	9	-84.64806	187.6114	10.525683	5.053890e-03
10	-6.578647	-0.3485384	4.31e-05	-0.0081510	-6.0e-07	-0.0083501	0.0002912	NA	0.0118731	NA	-0.0045846	NA	10	-87.51374	195.4134	18.327755	1.021945e-04
18	-6.619830	-0.2817268	4.22e-05	-0.0066989	-6.0e-07	-0.0111292	0.0029242	NA	0.0095350	NA	NA	-0.0038774	10	-87.72589	195.8377	18.752067	8.265882e-05
1	-6.793837	-0.0071012	3.79e-05	-0.0003277	-8.0e-07	-0.0104952	0.0000638	-7.00e-06	NA	NA	NA	NA	9	-89.52020	197.3556	20.269961	3.869749e-05
24	-6.791457	-0.0095030	3.41e-05	-0.0004282	-7.0e-07	-0.0123313	0.0034424	NA	NA	NA	0.0029834	-0.0047902	10	-88.51791	197.4218	20.336105	3.743862e-05
4	-6.792362	-0.0119234	3.40e-05	-0.0002822	-9.0e-07	-0.0105858	0.0000320	NA	NA	2.0e-07	NA	NA	9	-91.92854	202.1723	25.086646	3.481392e-06
3	-6.546620	-0.3881417	7.70e-05	-0.0097911	-5.0e-07	-0.0110606	0.0002367	-5.34e-05	0.0140699	NA	NA	NA	10	-92.07196	204.5299	27.444213	1.071060e-06
26	-6.595359	-0.3203397	4.24e-05	-0.0075827	-6.0e-07	-0.0095563	0.0022257	NA	0.0109002	NA	-0.0024481	-0.0027530	11	-91.56060	205.5852	28.499491	6.319213e-07
9	-6.795188	-0.0035460	4.36e-05	-0.0003527	-8.0e-07	-0.0114750	0.0000104	-1.58e-05	NA	NA	0.0015187	NA	10	-93.28777	206.9615	29.875824	3.175383e-07
17	-6.790340	-0.0127444	2.80e-05	-0.0003852	-8.0e-07	-0.0103599	0.0028331	8.50e-06	NA	NA	NA	-0.0037868	10	-93.54185	207.4697	30.383992	2.462912e-07
6	-6.545265	-0.4185973	4.46e-05	-0.0085371	-2.7e-06	-0.0111960	-0.0000539	NA	0.0130245	4.0e-06	NA	NA	10	-94.64428	209.6745	32.588849	8.178437e-08
20	-6.791009	-0.0105540	3.26e-05	-0.0003928	-9.0e-07	-0.0102326	0.0025179	NA	NA	2.0e-07	NA	-0.0033237	10	-96.07841	212.5428	35.457096	1.949119e-08
12	-6.792376	-0.0118581	3.40e-05	-0.0002830	-9.0e-07	-0.0105609	0.0000341	NA	NA	2.0e-07	-0.0000348	NA	10	-96.12841	212.6428	35.557100	1.854055e-08
11	-6.546790	-0.3878181	7.72e-05	-0.0097879	-5.0e-07	-0.0111019	0.0002344	-5.38e-05	0.0140635	NA	0.0000644	NA	11	-95.84961	214.1632	37.077509	8.669023e-09
19	-6.556102	-0.3736380	7.09e-05	-0.0093945	-5.0e-07	-0.0109759	0.0014473	-4.46e-05	0.0134426	NA	NA	-0.0016659	11	-96.24076	214.9455	37.859814	5.862662e-09
25	-6.793167	-0.0044397	4.12e-05	-0.0004673	-7.0e-07	-0.0129949	0.0033627	-1.18e-05	NA	NA	0.0041352	-0.0047097	11	-97.16621	216.7964	39.710712	2.323685e-09
14	-6.490797	-0.5026104	4.45e-05	-0.0105322	-2.7e-06	-0.0079572	0.0002189	NA	0.0159988	4.0e-06	-0.0047149	NA	11	-98.48958	219.4431	42.357454	6.186494e-10
22	-6.529743	-0.4396690	4.35e-05	-0.0091370	-2.8e-06	-0.0107908	0.0030334	NA	0.0137539	4.2e-06	NA	-0.0041341	11	-98.67549	219.8150	42.729279	5.136933e-10
5	-6.794334	-0.0057058	3.83e-05	-0.0003407	-7.0e-07	-0.0105065	0.0000685	-7.60e-06	NA	-1.0e-07	NA	NA	10	-100.92841	222.2428	45.157109	1.525834e-10
28	-6.788932	-0.0162012	3.38e-05	-0.0003694	-1.1e-06	-0.0124711	0.0035496	NA	NA	7.0e-07	0.0033702	-0.0049825	11	-99.96767	222.3993	45.313628	1.410976e-10
27	-6.561226	-0.3645445	7.41e-05	-0.0092555	-5.0e-07	-0.0117654	0.0016335	-4.98e-05	0.0131979	NA	0.0012566	-0.0019849	12	-99.92687	224.4030	47.317360	5.181015e-11
7	-6.489092	-0.4927366	7.38e-05	-0.0111805	-2.2e-06	-0.0108759	0.0001476	-4.68e-05	0.0166261	3.2e-06	NA	NA	11	-103.17758	228.8191	51.733441	5.694752e-12
30	-6.506666	-0.4760139	4.37e-05	-0.0099741	-2.8e-06	-0.0092833	0.0023628	NA	0.0150450	4.2e-06	-0.0023486	-0.0030542	12	-102.51344	229.5762	52.490494	3.900164e-12
13	-6.796232	-0.0006312	4.47e-05	-0.0003795	-6.0e-07	-0.0115476	0.0000170	-1.75e-05	NA	-2.0e-07	0.0015968	NA	11	-104.68194	231.8278	54.742163	1.265142e-12
21	-6.787454	-0.0205873	2.53e-05	-0.0003187	-1.1e-06	-0.0102908	0.0029741	1.26e-05	NA	7.0e-07	NA	-0.0040147	11	-104.90431	232.2726	55.186921	1.012888e-12
15	-6.483098	-0.5039427	6.92e-05	-0.0113082	-2.3e-06	-0.0100458	0.0001901	-3.95e-05	0.0168697	3.3e-06	-0.0012815	NA	12	-106.92194	238.3932	61.307488	4.747851e-14
23	-6.496504	-0.4833026	6.45e-05	-0.0107540	-2.4e-06	-0.0107296	0.0019254	-3.29e-05	0.0159927	3.5e-06	NA	-0.0024605	12	-107.28539	239.1201	62.034394	3.301045e-14
29	-6.790799	-0.0108871	3.88e-05	-0.0004124	-1.0e-06	-0.0128884	0.0034655	-8.10e-06	NA	5.0e-07	0.0040550	-0.0048745	12	-108.52941	241.6081	64.522435	9.514371e-15
31	-6.497205	-0.4820471	6.49e-05	-0.0107357	-2.4e-06	-0.0108121	0.0019433	-3.35e-05	0.0159601	3.5e-06	0.0001301	-0.0024913	13	-110.96642	248.5748	71.489128	2.921340e-16

Run the top model and 100 models for each one:

#In both cases the top model is 1/2/3/4/5/6 no interaction terms. With no models within delta < 4: 

#Run model for DR
Armenta.MCC.DR.top <- gls(MCC.DR ~ abs.Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = FALSE), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.DR.top, 'data/Armenta.MCC.DR.top.rds')

#Run model for ND
Armenta.MCC.ND.top <- gls(MCC.ND ~ abs.Colour.discriminability
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Armenta.tree, fixed = TRUE), #lambda = 1
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.ND.top, 'data/Armenta.MCC.ND.top.rds')

Armenta.data$abs.Colour.discriminability <- abs(Armenta.data$Colour.discriminability)

#Run the 100 models for DR and ND using the best model:
Armenta.data.noMCC <- inner_join(restricted.data %>% dplyr::select(TipLabel, binomial, SDi, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), Armenta.data, by = "binomial") %>% filter(abs.Colour.discriminability != "NA")
#Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
Armenta.DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, abs.Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "DR")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.DR.model.data <- lapply(Armenta.DR.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
Armenta.pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(MCC.Armenta.model.data$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
Armenta.DR.pgls.models <- mcmapply(function(x,y) {
  gls(DR ~ abs.Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    correlation = corPagel(1, phy = y, fixed = FALSE), 
    data = x, 
    method = "REML")
}, x = Armenta.DR.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.DR.pgls.models, "data/Armenta.DR.pgls.models.rds")

#Now for Node Density:
Armenta.ND.model.data <- lapply(nd.list, function(x) {
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, abs.Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.ND.model.data <- lapply(Armenta.ND.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Armenta.ND.pgls.models <- mcmapply(function(x,y) {
gls(ND ~ abs.Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(1, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = Armenta.ND.model.data, y = Armenta.pruned.trees, 
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.ND.pgls.models, "data/Armenta.ND.pgls.models.rds")

#Run the BAMM models

Armenta.MCC.Lambda.top <- gls(log(mean.lambda) ~ abs.Colour.discriminability
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ I(1/log(var.lambda)),
                correlation = corBrownian(phy = pruned.MCC.Armenta.tree), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.Lambda.top, 'data/Armenta.MCC.Lambda.top.rds')

Armenta.MCC.Mu.top <- gls(log(mean.mu) ~ abs.Colour.discriminability
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ I(1/log(var.mu)),
                correlation = corBrownian(phy = pruned.MCC.Armenta.tree), 
                data = MCC.Armenta.model.data, 
                method = "REML")
saveRDS(Armenta.MCC.Mu.top, 'data/Armenta.MCC.Mu.top.rds')

Armenta.BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Armenta.data.noMCC %>% dplyr::select(binomial, TipLabel, SDi, abs.Colour.discriminability, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Armenta.BAMM.model.data <- lapply(Armenta.BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Armenta.BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(log(mean.lambda) ~ abs.Colour.discriminability 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ I(1/log(var.lambda)),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Armenta.BAMM.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.BAMM.lambda.pgls.models, "data/Armenta.BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
Armenta.BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(log(mean.mu) ~ abs.Colour.discriminability
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ I(1/log(var.mu)),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Armenta.BAMM.model.data, y = Armenta.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Armenta.BAMM.mu.pgls.models, "data/Armenta.BAMM.mu.pgls.models.rds")

Using the MCC models as well as the 100 PGLS models we can generate the plots similar to Supp Fig S8:

Armenta.DR.pgls.models <- readRDS("data/Armenta.DR.pgls.models.rds")
Armenta.MCC.DR.top <- readRDS('data/Armenta.MCC.DR.top.rds')

Armenta.DR.pgls.summary <- lapply(Armenta.DR.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.DR.summary <- data.frame(Armenta.MCC.DR.top$coefficients, confint(Armenta.MCC.DR.top)) %>% tibble::rownames_to_column()

Armenta.DR.pgls.summary <- bind_rows(Armenta.DR.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.DR.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `abs.Colour.discriminability` = "Sexual Dichromatism"
                    )

Armenta.DR.pgls.summary$Parameter = factor(Armenta.DR.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'abs.Colour.discriminability'))

Armenta.MCC.DR.summary$Parameter = factor(Armenta.MCC.DR.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'abs.Colour.discriminability'))

Armenta.DR.plot <-Armenta.DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.DR.plot <- Armenta.DR.plot + geom_errorbarh(data = Armenta.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta DR")

####################################
#For ND

Armenta.ND.pgls.models <- readRDS("data/Armenta.ND.pgls.models.rds")
Armenta.MCC.ND.top <- readRDS('data/Armenta.MCC.ND.top.rds')

Armenta.ND.pgls.summary <- lapply(Armenta.ND.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.ND.summary <- data.frame(Armenta.MCC.ND.top$coefficients, confint(Armenta.MCC.ND.top)) %>% tibble::rownames_to_column()

Armenta.ND.pgls.summary <- bind_rows(Armenta.ND.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.ND.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.ND.pgls.summary$Parameter = factor(Armenta.ND.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'abs.Colour.discriminability'))

Armenta.MCC.ND.summary$Parameter = factor(Armenta.MCC.ND.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'abs.Colour.discriminability'))

Armenta.ND.plot <-Armenta.ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.ND.plot <- Armenta.ND.plot + geom_errorbarh(data = Armenta.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta ND")

######################

#For Lambda
Armenta.BAMM.lambda.pgls.models <- readRDS("data/Armenta.BAMM.lambda.pgls.models.rds")
Armenta.MCC.Lambda.top <- readRDS('data/Armenta.MCC.Lambda.top.rds')

Armenta.Lambda.pgls.summary <- lapply(Armenta.BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.Lambda.summary <- data.frame(Armenta.MCC.Lambda.top$coefficients, confint(Armenta.MCC.Lambda.top)) %>% tibble::rownames_to_column()

Armenta.Lambda.pgls.summary <- bind_rows(Armenta.Lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.Lambda.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.Lambda.pgls.summary$Parameter = factor(Armenta.Lambda.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'abs.Colour.discriminability'))

Armenta.MCC.Lambda.summary$Parameter = factor(Armenta.MCC.Lambda.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'abs.Colour.discriminability'))

Armenta.Lambda.plot <-Armenta.Lambda.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.Lambda.plot <- Armenta.Lambda.plot + geom_errorbarh(data = Armenta.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta BAMM Speciaton")

######################

#For Mu

Armenta.BAMM.Mu.pgls.models <- readRDS("data/Armenta.BAMM.mu.pgls.models.rds")
Armenta.MCC.Mu.top <- readRDS('data/Armenta.MCC.Mu.top.rds')

Armenta.Mu.pgls.summary <- lapply(Armenta.BAMM.Mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Armenta.MCC.Mu.summary <- data.frame(Armenta.MCC.Mu.top$coefficients, confint(Armenta.MCC.Mu.top)) %>% tibble::rownames_to_column()

Armenta.Mu.pgls.summary <- bind_rows(Armenta.Mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Armenta.MCC.Mu.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Armenta.Mu.pgls.summary$Parameter = factor(Armenta.Mu.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'abs.Colour.discriminability'))

Armenta.MCC.Mu.summary$Parameter = factor(Armenta.MCC.Mu.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'abs.Colour.discriminability'))

Armenta.Mu.plot <-Armenta.Mu.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Armenta.Mu.plot <- Armenta.Mu.plot + geom_errorbarh(data = Armenta.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Armenta.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Armenta BAMM Extinction")

#########################################

Plot the results

grid.arrange(symmetrise_scale(Armenta.DR.plot, "x"),
             symmetrise_scale(Armenta.ND.plot, "x"),
             symmetrise_scale(Armenta.Lambda.plot, "x"), 
             symmetrise_scale(Armenta.Mu.plot, "x"), 
             nrow = 1)

Figure S9: Measures of sexual dichromatism using a more precise measure in spectrophotometry do not change the main patterns seen within the PGLS models analysed here. This dataset is a subset of the complete dataset (n = 581), thus drawing conclusions for the other predictors (e.g. borderline long term temperature variation and spatial temperature variation) potentially risks type I error; additionally these environmental predictors were not measured differently in this analysis compared to the previous above so we have no obvious reason to make conclusions from this subset for the effects of the environmental predictors.

Table S9: Highest posterior density intervals are calculated for the above figure (ie; models using Spectrophotometry values of sexual dichromatism). The HPD range is determined using the hdi function with a 95 % credible interval. These intervals do not take into account the variance associated with each interval and thus are not an estimate of model precision. Intervals not overlapping zero suggest that 95 % of trees from the posterior generate a model estimate for the given parameter that in the same direction (+ or -). These intervals are calculated in the same way as in Table S7.

Armenta.hpd.DR.top <- list()
Armenta.DR.pgls.summary <- na.omit(Armenta.DR.pgls.summary)
for(x in unique(Armenta.DR.pgls.summary$Parameter)) {
Armenta.hpd.DR.top[[x]] = hdi(Armenta.DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.DR.top <- bind_rows(Armenta.hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper"))


saveRDS(Armenta.hpd.DR.top, 'data/Armenta.hpd.DR.top.rds')

#For ND
Armenta.hpd.ND.top <- list()
Armenta.ND.pgls.summary <- na.omit(Armenta.ND.pgls.summary)
for(x in unique(Armenta.ND.pgls.summary$Parameter)) {
Armenta.hpd.ND.top[[x]] = hdi(Armenta.ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.ND.top <- bind_rows(Armenta.hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Armenta.hpd.ND.top, 'data/Armenta.hpd.ND.top.rds')

###############
Armenta.hpd.Lambda.top <- list()
Armenta.Lambda.pgls.summary <- na.omit(Armenta.Lambda.pgls.summary)
for(x in unique(Armenta.Lambda.pgls.summary$Parameter)) {
Armenta.hpd.Lambda.top[[x]] = hdi(Armenta.Lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.Lambda.top <- bind_rows(Armenta.hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) 
saveRDS(Armenta.hpd.Lambda.top, 'data/Armenta.hpd.Lambda.top.rds')

################
Armenta.hpd.Mu.top <- list()
Armenta.Mu.pgls.summary <- na.omit(Armenta.Mu.pgls.summary)
for(x in unique(Armenta.Mu.pgls.summary$Parameter)) {
Armenta.hpd.Mu.top[[x]] = hdi(Armenta.Mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Armenta.hpd.Mu.top <- bind_rows(Armenta.hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Armenta.hpd.Mu.top, 'data/Armenta.hpd.Mu.top.rds')

Armenta.hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta DR HPD Interval")

Armenta DR HPD Interval
	abs.Colour.discriminability	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.0393	-0.0143	-0.000139	0.0108	-0.0427	-4.27e-06
Upper	0.0421	0.0316	0.000159	0.0347	0.0076	1.44e-05

Armenta.hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta ND HPD Interval")

Armenta ND HPD Interval
	abs.Colour.discriminability	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.0021	-0.000828	-4.96e-06	0.000491	-0.00228	-6.59e-08
Upper	0.0026	0.00121	1.38e-05	0.00181	0.000572	6.88e-07

Armenta.hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta BAMM Speciation HPD Interval")

Armenta BAMM Speciation HPD Interval
	abs.Colour.discriminability	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.0243	-0.00942	-6.03e-05	-0.00231	-0.0114	-3.9e-06
Upper	0.0144	0.00725	2.65e-05	0.00543	0.00317	3.2e-06

Armenta.hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "Armenta BAMM Extinction HPD Interval")

Armenta BAMM Extinction HPD Interval
	abs.Colour.discriminability	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.0337	-0.0114	-5.71e-05	-0.00565	-0.0228	-5.39e-06
Upper	0.0269	0.0112	7.6e-05	0.00487	0.00156	4.43e-06

Analysis using male-bias measure of sexual selection

Sexual dichromatism is expected to be a good measure of sexual selection strength in birds. However the relationship between sexual dichromatism and male-biased measures of sexual selection (social mating system, sexual size dimorphism, and paternal care) is expected to be relatively noisy given the precision of the measurment used (Dale et al. 2015). Here we use a dataset of sexual selection for 2,465 species within the orginal dataset of sexual dichromatism scores from RGB measures. This male-bias sexual selection score is based on three components, combined in a phylogenetic pca (ppca) with the following loadings:

sexual size dimorphism: 0.37
social polygyny: 0.57
paternal care: -0.57

As such, high values for this score are expected to have high sexual selection (high dimorphism, high polygyny and low paternal care). PGLS models using this dataset are calculated in the same process as above.

SS.subset <- plumage.scores %>% filter(Sexual_selection_ppca != "NA")
SS.subset <- inner_join(SS.subset, restricted.data %>% dplyr::select(TipLabel, binomial, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP, MCC.DR, MCC.ND), by = "TipLabel")
SS.subset <- inner_join(SS.subset, MCC.BAMM.df, by = "TipLabel")

SS.subset %>% ggplot(aes(x = SDi, y = Sexual_selection_ppca))+
  geom_point()+
  geom_smooth(method = "lm")+
  theme_minimal()

Figure S10: The relationship between absolute sexual dichromatism and the sexual selection. Although positively correlated the distribution is very noisy with species having high sexual dichromatism and low sexual selection and visa versa.

#Prune tree
pruned.MCC.Subset.tree <- drop.tip(MCC.passerine,MCC.passerine$tip.label[-match(SS.subset$TipLabel, MCC.passerine$tip.label)])

#Set rownames to match tree
rownames(SS.subset) <- SS.subset$TipLabel

#Run a corPagel model to estimate lambda for DR
MCC.DR.Subset.global <- gls(MCC.DR ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = FALSE), 
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.ND.Subset.global <- gls(MCC.ND ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), #lambda = 1
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Lambda.Subset.global <- gls(log(mean.lambda) ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                weights = ~ I(1/log(var.lambda)),
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), 
                data = SS.subset, 
                method = "REML")

#Run a corPagel model to estimate lambda for DR
MCC.Extinction.Subset.global <- gls(log(mean.mu) ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP
                         + Sexual_selection_ppca*log(range.size.m2)
                         + Sexual_selection_ppca*bioclim4
                         + Sexual_selection_ppca*residuals.PC1
                         + Sexual_selection_ppca*PC1.LIG
                         + Sexual_selection_ppca*NPP,
                weights = ~ I(1/log(var.mu)),
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), 
                data = SS.subset, 
                method = "REML")

#Set up cluster
cores<-8
clusterType <- if(length(find.package("snow", quiet = TRUE))) "SOCK" else "PSOCK"
clust <- try(makeCluster(getOption("cl.cores", cores), type = clusterType))
myclust<-clust

#Export data and packages to cluster
clusterExport(myclust, c("SS.subset"), envir=environment())
clusterExport(myclust, c("pruned.MCC.Subset.tree"), envir=environment())
clusterEvalQ(myclust, library(nlme))
clusterEvalQ(myclust, library(ape))
clusterEvalQ(myclust, library(MuMIn))

#Dredged models:

Subset.dredged.ND.model <- pdredge(MCC.ND.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.ND.model, 'data/Subset.dredged.ND.model.rds')

Subset.dredged.DR.model <- pdredge(MCC.DR.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.DR.model, 'data/Subset.dredged.DR.model.rds')

Subset.dredged.spec.model <- pdredge(MCC.Lambda.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.spec.model, 'data/Subset.dredged.spec.model.rds')

Subset.dredged.extinct.model <- pdredge(MCC.Extinction.Subset.global, fixed = c("Sexual_selection_ppca", "log(range.size.m2)", "bioclim4", "residuals.PC1", "PC1.LIG", "NPP"), cluster=myclust)

saveRDS(Subset.dredged.spec.model, 'data/Subset.dredged.extinct.model.rds')

Table S10: All top models (\(\delta AICc > 4\)) using male-biased sexual selection measures do not contain interactions between sexual selection and the environmental variables. Thus interaction terms are not included in further analysis.

Subset.dredged.ND.model <- readRDS("data/Subset.dredged.ND.model.rds")
Subset.dredged.DR.model <- readRDS("data/Subset.dredged.DR.model.rds")
Subset.dredged.spec.model <- readRDS("data/Subset.dredged.spec.model.rds")
Subset.dredged.extinct.model <- readRDS("data/Subset.dredged.extinct.model.rds")


kable(Subset.dredged.ND.model, "html", caption = "ND Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

ND Male-bias Sexual Selection Dredge Table
	(Intercept)	bioclim4	log(range.size.m2)	PC1.LIG	residuals.PC1	Sexual_selection_ppca	bioclim4:Sexual_selection_ppca	log(range.size.m2):Sexual_selection_ppca	NPP:Sexual_selection_ppca	PC1.LIG:Sexual_selection_ppca	residuals.PC1:Sexual_selection_ppca	df	logLik	AICc	delta	weight
0	0.1484807	-5e-07	-0.0001800	0.0001073	0.0000788	0.0004381	NA	NA	NA	NA	NA	8	5451.548	-10887.04	0.00000	9.994587e-01
16	0.1482246	-5e-07	-0.0001720	0.0000929	0.0000644	0.0003782	NA	NA	NA	NA	-0.0001167	9	5444.064	-10870.05	16.98341	2.050523e-04
8	0.1479484	-5e-07	-0.0001570	0.0000777	0.0001046	0.0004928	NA	NA	NA	0.0001280	NA	9	5444.026	-10869.98	17.05889	1.974574e-04
2	0.1483243	-5e-07	-0.0001765	0.0001017	0.0000834	0.0018018	NA	-4.87e-05	NA	NA	NA	9	5443.645	-10869.22	17.82105	1.348878e-04
1	0.1478764	-9e-07	-0.0001561	0.0000738	0.0000871	0.0012068	-2.3e-06	NA	NA	NA	NA	9	5440.065	-10862.06	24.98016	3.761776e-06
4	0.1484220	-6e-07	-0.0001730	0.0001024	0.0000814	0.0008188	NA	NA	0e+00	NA	NA	9	5436.296	-10854.52	32.51872	8.677880e-08
24	0.1477981	-5e-07	-0.0001528	0.0000692	0.0000925	0.0004458	NA	NA	NA	0.0001196	-0.0000845	10	5436.455	-10852.82	34.21679	3.712640e-08
18	0.1479767	-6e-07	-0.0001662	0.0000832	0.0000688	0.0022340	NA	-6.65e-05	NA	NA	-0.0001323	10	5436.226	-10852.36	34.67399	2.953948e-08
10	0.1478037	-5e-07	-0.0001538	0.0000725	0.0001088	0.0017783	NA	-4.59e-05	NA	0.0001274	NA	10	5436.116	-10852.14	34.89512	2.644759e-08
17	0.1478092	-9e-07	-0.0001543	0.0000700	0.0000806	0.0011337	-2.2e-06	NA	NA	NA	-0.0000480	10	5432.459	-10844.83	42.20969	6.824126e-10
9	0.1474813	-8e-07	-0.0001389	0.0000518	0.0001085	0.0011728	-2.1e-06	NA	NA	0.0001104	NA	10	5432.373	-10844.66	42.38150	6.262388e-10
3	0.1478909	-9e-07	-0.0001563	0.0000742	0.0000866	0.0010291	-2.3e-06	6.70e-06	NA	NA	NA	10	5432.150	-10844.21	42.82661	5.012840e-10
20	0.1482076	-6e-07	-0.0001673	0.0000906	0.0000682	0.0006967	NA	NA	0e+00	NA	-0.0001026	10	5428.776	-10837.46	49.57473	1.716915e-11
12	0.1478969	-6e-07	-0.0001505	0.0000733	0.0001068	0.0008569	NA	NA	0e+00	0.0001269	NA	10	5428.762	-10837.43	49.60299	1.692825e-11
6	0.1482266	-6e-07	-0.0001680	0.0000951	0.0000872	0.0025093	NA	-5.90e-05	0e+00	NA	NA	10	5428.428	-10836.77	50.27032	1.212563e-11
26	0.1475846	-6e-07	-0.0001479	0.0000610	0.0000958	0.0021036	NA	-5.95e-05	NA	0.0001173	-0.0000990	11	5428.596	-10835.09	51.95195	5.230492e-12
5	0.1477933	-1e-06	-0.0001475	0.0000674	0.0000903	0.0016611	-2.4e-06	NA	-1e-07	NA	NA	10	5424.857	-10829.62	57.41215	3.410950e-13
25	0.1474541	-8e-07	-0.0001383	0.0000503	0.0001049	0.0011368	-2.0e-06	NA	NA	0.0001086	-0.0000240	11	5424.751	-10827.39	59.64290	1.118085e-13
19	0.1477962	-9e-07	-0.0001541	0.0000696	0.0000808	0.0012568	-2.1e-06	-4.70e-06	NA	NA	-0.0000499	11	5424.591	-10827.07	59.96353	9.524668e-14
11	0.1474880	-8e-07	-0.0001390	0.0000520	0.0001082	0.0010942	-2.1e-06	2.90e-06	NA	0.0001103	NA	11	5424.457	-10826.81	60.23059	8.334079e-14
28	0.1477791	-6e-07	-0.0001479	0.0000668	0.0000965	0.0007717	NA	NA	0e+00	0.0001201	-0.0000699	11	5421.172	-10820.24	66.80118	3.119414e-15
22	0.1479328	-7e-07	-0.0001603	0.0000797	0.0000735	0.0027753	NA	-7.32e-05	0e+00	NA	-0.0001181	11	5420.966	-10819.82	67.21289	2.539046e-15
14	0.1477158	-6e-07	-0.0001460	0.0000665	0.0001122	0.0024557	NA	-5.58e-05	0e+00	0.0001260	NA	11	5420.885	-10819.66	67.37476	2.341650e-15
21	0.1477617	-1e-06	-0.0001470	0.0000656	0.0000868	0.0016033	-2.3e-06	NA	0e+00	NA	-0.0000251	11	5417.249	-10812.39	74.64696	6.171400e-17
13	0.1474084	-1e-06	-0.0001309	0.0000460	0.0001112	0.0016076	-2.1e-06	NA	0e+00	0.0001085	NA	11	5417.149	-10812.19	74.84702	5.583943e-17
7	0.1477851	-1e-06	-0.0001473	0.0000671	0.0000906	0.0017584	-2.4e-06	-3.60e-06	-1e-07	NA	NA	11	5416.948	-10811.79	75.24841	4.568560e-17
27	0.1474468	-8e-07	-0.0001381	0.0000501	0.0001050	0.0012077	-2.0e-06	-2.70e-06	NA	0.0001085	-0.0000251	12	5416.883	-10809.64	77.39880	1.558939e-17
30	0.1475395	-6e-07	-0.0001419	0.0000575	0.0001007	0.0026498	NA	-6.62e-05	0e+00	0.0001176	-0.0000846	12	5413.338	-10802.55	84.48797	4.502300e-19
29	0.1474072	-1e-06	-0.0001309	0.0000459	0.0001110	0.0016049	-2.1e-06	NA	0e+00	0.0001084	-0.0000012	12	5409.540	-10794.95	92.08377	1.009318e-20
23	0.1477348	-1e-06	-0.0001465	0.0000647	0.0000871	0.0018596	-2.2e-06	-9.70e-06	0e+00	NA	-0.0000288	12	5409.384	-10794.64	92.39589	8.634816e-21
15	0.1473924	-1e-06	-0.0001306	0.0000455	0.0001118	0.0017935	-2.1e-06	-6.80e-06	0e+00	0.0001087	NA	12	5409.241	-10794.35	92.68281	7.480788e-21
31	0.1473862	-1e-06	-0.0001305	0.0000452	0.0001112	0.0018080	-2.1e-06	-7.70e-06	0e+00	0.0001084	-0.0000042	13	5401.675	-10777.20	109.83562	1.410139e-24

kable(Subset.dredged.DR.model, "html", caption = "DR Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

DR Male-bias Sexual Selection Dredge Table
	(Intercept)	bioclim4	log(range.size.m2)	NPP	PC1.LIG	residuals.PC1	Sexual_selection_ppca	bioclim4:Sexual_selection_ppca	log(range.size.m2):Sexual_selection_ppca	NPP:Sexual_selection_ppca	PC1.LIG:Sexual_selection_ppca	residuals.PC1:Sexual_selection_ppca	df	logLik	AICc	delta	weight
0	-2.644369	2.45e-05	-0.0113055	-7.0e-07	0.0068813	0.0004409	0.0388722	NA	NA	NA	NA	NA	9	-1717.437	3452.947	0.00000	9.851357e-01
16	-2.639661	2.53e-05	-0.0115625	-7.0e-07	0.0066554	-0.0000035	0.0381221	NA	NA	NA	NA	-0.0043664	10	-1721.450	3462.989	10.04197	6.499963e-03
2	-2.651446	2.27e-05	-0.0111253	-7.0e-07	0.0068731	0.0003113	0.1032726	NA	-0.0023191	NA	NA	NA	10	-1721.901	3463.892	10.94493	4.138415e-03
8	-2.643685	2.47e-05	-0.0113672	-8.0e-07	0.0067933	0.0002826	0.0391218	NA	NA	NA	-0.0017643	NA	10	-1721.915	3463.919	10.97215	4.082471e-03
1	-2.644777	1.73e-05	-0.0112474	-8.0e-07	0.0069016	0.0002884	0.0537640	-4.40e-05	NA	NA	NA	NA	10	-1726.048	3472.186	19.23901	6.543319e-05
18	-2.648357	2.30e-05	-0.0113613	-7.0e-07	0.0066118	-0.0002393	0.1227636	NA	-0.0030524	NA	NA	-0.0050116	11	-1725.773	3473.654	20.70712	3.140527e-05
24	-2.639316	2.54e-05	-0.0115972	-7.0e-07	0.0066022	-0.0000984	0.0383393	NA	NA	NA	-0.0012567	-0.0041819	11	-1725.963	3474.034	21.08668	2.597662e-05
10	-2.650479	2.30e-05	-0.0111903	-7.0e-07	0.0067936	0.0001735	0.1006675	NA	-0.0022170	NA	-0.0015983	NA	11	-1726.392	3474.892	21.94541	1.690876e-05
4	-2.639136	2.37e-05	-0.0113014	-1.1e-06	0.0067364	0.0005491	0.0550570	NA	NA	-1.6e-06	NA	NA	10	-1729.238	3478.565	25.61792	2.695472e-06
17	-2.641237	1.95e-05	-0.0114476	-8.0e-07	0.0067343	0.0000040	0.0497120	-3.36e-05	NA	NA	NA	-0.0031418	11	-1730.243	3482.593	29.64645	3.596263e-07
9	-2.644302	1.77e-05	-0.0112841	-8.0e-07	0.0068566	0.0002153	0.0532073	-4.20e-05	NA	NA	-0.0008810	NA	11	-1730.569	3483.246	30.29926	2.594752e-07
3	-2.648526	1.70e-05	-0.0111570	-7.0e-07	0.0068950	0.0002315	0.0876844	-3.97e-05	-0.0012746	NA	NA	NA	11	-1730.598	3483.303	30.35634	2.521738e-07
26	-2.647793	2.32e-05	-0.0113950	-7.0e-07	0.0065727	-0.0003052	0.1206409	NA	-0.0029697	NA	-0.0009490	-0.0048550	12	-1730.301	3484.728	31.78144	1.236643e-07
20	-2.635570	2.45e-05	-0.0115406	-1.0e-06	0.0065504	0.0001210	0.0518306	NA	NA	-1.4e-06	NA	-0.0040359	11	-1733.315	3488.737	35.78960	1.666795e-08
6	-2.646854	2.14e-05	-0.0110830	-1.1e-06	0.0067038	0.0004113	0.1349625	NA	-0.0027849	-1.9e-06	NA	NA	11	-1733.615	3489.338	36.39048	1.234251e-08
12	-2.638336	2.39e-05	-0.0113656	-1.1e-06	0.0066414	0.0003852	0.0556198	NA	NA	-1.7e-06	-0.0018515	NA	11	-1733.706	3489.519	36.57214	1.127085e-08
19	-2.647224	1.97e-05	-0.0113350	-7.0e-07	0.0066756	-0.0001778	0.1100270	-2.27e-05	-0.0023104	NA	NA	-0.0040281	12	-1734.642	3493.410	40.46326	1.610686e-09
25	-2.640864	1.99e-05	-0.0114754	-8.0e-07	0.0067004	-0.0000514	0.0493311	-3.21e-05	NA	NA	-0.0007139	-0.0030923	12	-1734.769	3493.665	40.71787	1.418151e-09
11	-2.648042	1.75e-05	-0.0111939	-8.0e-07	0.0068500	0.0001586	0.0870861	-3.76e-05	-0.0012730	NA	-0.0008797	NA	12	-1735.118	3494.364	41.41693	9.998225e-10
5	-2.639640	1.65e-05	-0.0112442	-1.1e-06	0.0067596	0.0003957	0.0694708	-4.36e-05	NA	-1.6e-06	NA	NA	11	-1737.859	3497.825	44.87789	1.771687e-10
22	-2.644516	2.18e-05	-0.0113100	-1.1e-06	0.0064809	-0.0001171	0.1492103	NA	-0.0034126	-1.7e-06	NA	-0.0046899	12	-1737.559	3499.245	46.29775	8.710986e-11
28	-2.635111	2.47e-05	-0.0115772	-1.1e-06	0.0064896	0.0000208	0.0524230	NA	NA	-1.4e-06	-0.0013747	-0.0038249	12	-1737.818	3499.763	46.81566	6.723639e-11
14	-2.645817	2.17e-05	-0.0111498	-1.2e-06	0.0066197	0.0002691	0.1325514	NA	-0.0026832	-1.9e-06	-0.0016616	NA	12	-1738.100	3500.327	47.37978	5.071161e-11
27	-2.646817	2.00e-05	-0.0113622	-7.0e-07	0.0066447	-0.0002281	0.1093216	-2.14e-05	-0.0022967	NA	-0.0006587	-0.0039771	13	-1739.168	3504.485	51.53820	6.340417e-12
21	-2.637030	1.86e-05	-0.0114218	-1.1e-06	0.0066273	0.0001336	0.0642539	-3.45e-05	NA	-1.4e-06	NA	-0.0027665	12	-1742.094	3508.315	55.36799	9.343053e-13
7	-2.644376	1.61e-05	-0.0111162	-1.1e-06	0.0067357	0.0003290	0.1181574	-3.75e-05	-0.0017672	-1.8e-06	NA	NA	12	-1742.351	3508.828	55.88132	7.228072e-13
13	-2.639083	1.70e-05	-0.0112841	-1.1e-06	0.0067082	0.0003160	0.0690252	-4.13e-05	NA	-1.6e-06	-0.0009791	NA	12	-1742.374	3508.875	55.92791	7.061633e-13
30	-2.643870	2.20e-05	-0.0113459	-1.1e-06	0.0064361	-0.0001880	0.1471758	NA	-0.0033252	-1.7e-06	-0.0010505	-0.0045120	13	-1742.080	3510.308	57.36050	3.450013e-13
23	-2.643470	1.87e-05	-0.0112857	-1.1e-06	0.0065443	-0.0000596	0.1365454	-2.19e-05	-0.0026920	-1.6e-06	NA	-0.0037463	13	-1746.434	3519.016	66.06877	4.434500e-15
29	-2.636582	1.90e-05	-0.0114524	-1.1e-06	0.0065870	0.0000719	0.0640021	-3.28e-05	NA	-1.4e-06	-0.0008236	-0.0027038	13	-1746.614	3519.376	66.42937	3.702894e-15
15	-2.643820	1.66e-05	-0.0111564	-1.1e-06	0.0066838	0.0002487	0.1177972	-3.52e-05	-0.0017704	-1.8e-06	-0.0009854	NA	13	-1746.865	3519.878	66.93090	2.881620e-15
31	-2.642995	1.90e-05	-0.0113159	-1.1e-06	0.0065069	-0.0001168	0.1359737	-2.04e-05	-0.0026794	-1.6e-06	-0.0007727	-0.0036828	14	-1750.955	3530.081	77.13422	1.753934e-17

kable(Subset.dredged.spec.model, "html", caption = "BAMM Speciation Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

BAMM Speciation Male-bias Sexual Selection Dredge Table
	(Intercept)	bioclim4	log(range.size.m2)	PC1.LIG	residuals.PC1	Sexual_selection_ppca	bioclim4:Sexual_selection_ppca	log(range.size.m2):Sexual_selection_ppca	NPP:Sexual_selection_ppca	PC1.LIG:Sexual_selection_ppca	residuals.PC1:Sexual_selection_ppca	df	logLik	AICc	delta	weight
0	-2.835512	-7.3e-06	0.0004994	-0.0003236	0.0004112	-0.0003583	NA	NA	NA	NA	NA	8	2617.018	-5217.977	0.00000	9.982022e-01
8	-2.837029	-7.4e-06	0.0005575	-0.0003517	0.0004687	-0.0004496	NA	NA	NA	0.0004539	NA	9	2610.740	-5203.407	14.56967	6.846268e-04
16	-2.835892	-7.3e-06	0.0005110	-0.0003416	0.0003737	-0.0004570	NA	NA	NA	NA	-0.0002267	9	2610.652	-5203.231	14.74603	6.268432e-04
2	-2.835374	-7.2e-06	0.0004942	-0.0003057	0.0003976	-0.0048344	NA	0.0001612	NA	NA	NA	9	2610.384	-5202.694	15.28266	4.793257e-04
1	-2.835131	-7.0e-06	0.0004827	-0.0003044	0.0004030	-0.0008253	1.5e-06	NA	NA	NA	NA	9	2605.952	-5193.830	24.14720	5.697974e-06
24	-2.837371	-7.4e-06	0.0005678	-0.0003683	0.0004333	-0.0005410	NA	NA	NA	0.0004501	-0.0002114	10	2604.367	-5188.644	29.33283	4.262599e-07
10	-2.836889	-7.3e-06	0.0005522	-0.0003347	0.0004556	-0.0046804	NA	0.0001523	NA	0.0004513	NA	10	2604.101	-5188.112	29.86480	3.267087e-07
18	-2.835724	-7.3e-06	0.0005051	-0.0003243	0.0003679	-0.0041362	NA	0.0001330	NA	NA	-0.0001943	10	2604.018	-5187.947	30.03002	3.008033e-07
4	-2.835551	-7.3e-06	0.0005019	-0.0003257	0.0004112	-0.0001576	NA	NA	0	NA	NA	9	2602.884	-5187.695	30.28165	2.652418e-07
17	-2.835405	-6.9e-06	0.0004876	-0.0003165	0.0003449	-0.0012820	2.5e-06	NA	NA	NA	-0.0003181	10	2599.694	-5179.298	38.67931	3.982116e-09
9	-2.836626	-7.1e-06	0.0005397	-0.0003311	0.0004602	-0.0009541	1.6e-06	NA	NA	0.0004565	NA	10	2599.679	-5179.269	38.70807	3.925262e-09
3	-2.835149	-7.1e-06	0.0004842	-0.0002961	0.0003945	-0.0044284	1.0e-06	0.0001357	NA	NA	NA	10	2599.353	-5178.616	39.36092	2.832071e-09
26	-2.837205	-7.4e-06	0.0005620	-0.0003517	0.0004276	-0.0040318	NA	0.0001262	NA	0.0004485	-0.0001808	11	2597.730	-5173.352	44.62511	2.037042e-10
12	-2.837067	-7.4e-06	0.0005600	-0.0003538	0.0004688	-0.0002498	NA	NA	0	0.0004539	NA	10	2596.607	-5173.124	44.85298	1.817677e-10
20	-2.835901	-7.3e-06	0.0005117	-0.0003421	0.0003741	-0.0003955	NA	NA	0	NA	-0.0002247	10	2596.527	-5172.964	45.01281	1.678074e-10
6	-2.835390	-7.3e-06	0.0004952	-0.0003067	0.0003978	-0.0047065	NA	0.0001593	0	NA	NA	10	2596.258	-5172.427	45.55036	1.282580e-10
25	-2.836876	-6.9e-06	0.0005439	-0.0003425	0.0004039	-0.0013918	2.6e-06	NA	NA	0.0004525	-0.0003056	11	2593.413	-5164.719	53.25808	2.718774e-12
19	-2.835390	-6.9e-06	0.0004878	-0.0003117	0.0003459	-0.0028671	2.2e-06	0.0000611	NA	NA	-0.0002915	11	2593.127	-5164.146	53.83075	2.041833e-12
11	-2.836632	-7.1e-06	0.0005408	-0.0003235	0.0004522	-0.0042024	1.1e-06	0.0001223	NA	0.0004536	NA	11	2593.075	-5164.042	53.93520	1.937935e-12
5	-2.835167	-7.0e-06	0.0004849	-0.0003063	0.0004032	-0.0006555	1.5e-06	NA	0	NA	NA	10	2591.819	-5163.548	54.42891	1.514015e-12
28	-2.837380	-7.4e-06	0.0005686	-0.0003688	0.0004337	-0.0004705	NA	NA	0	0.0004502	-0.0002091	11	2590.242	-5158.376	59.60107	1.140276e-13
14	-2.836906	-7.3e-06	0.0005533	-0.0003357	0.0004558	-0.0045421	NA	0.0001503	0	0.0004513	NA	11	2589.975	-5157.843	60.13398	8.735578e-14
22	-2.835720	-7.3e-06	0.0005048	-0.0003240	0.0003677	-0.0041769	NA	0.0001335	0	NA	-0.0001950	11	2589.900	-5157.692	60.28495	8.100421e-14
27	-2.836861	-7.0e-06	0.0005441	-0.0003386	0.0004046	-0.0026825	2.3e-06	0.0000498	NA	0.0004516	-0.0002839	12	2586.845	-5149.562	68.41474	1.390423e-15
21	-2.835393	-6.9e-06	0.0004868	-0.0003158	0.0003443	-0.0013506	2.5e-06	NA	0	NA	-0.0003208	11	2585.575	-5149.042	68.93502	1.071937e-15
13	-2.836661	-7.1e-06	0.0005419	-0.0003330	0.0004604	-0.0007883	1.6e-06	NA	0	0.0004564	NA	11	2585.547	-5148.986	68.99156	1.042059e-15
7	-2.835165	-7.1e-06	0.0004852	-0.0002970	0.0003947	-0.0043083	1.0e-06	0.0001339	0	NA	NA	11	2585.227	-5148.347	69.62984	7.573414e-16
30	-2.837203	-7.4e-06	0.0005618	-0.0003516	0.0004275	-0.0040513	NA	0.0001264	0	0.0004485	-0.0001811	12	2583.611	-5143.095	74.88169	5.481093e-17
29	-2.836865	-6.9e-06	0.0005432	-0.0003419	0.0004034	-0.0014547	2.6e-06	NA	0	0.0004525	-0.0003080	12	2579.294	-5134.462	83.51548	7.312432e-19
23	-2.835374	-6.9e-06	0.0004867	-0.0003108	0.0003452	-0.0029967	2.2e-06	0.0000624	0	NA	-0.0002946	12	2579.012	-5133.896	84.08130	5.510563e-19
15	-2.836649	-7.1e-06	0.0005418	-0.0003245	0.0004524	-0.0040732	1.1e-06	0.0001205	0	0.0004536	NA	12	2578.949	-5133.771	84.20563	5.178431e-19
31	-2.836846	-7.0e-06	0.0005430	-0.0003377	0.0004039	-0.0027978	2.3e-06	0.0000509	0	0.0004516	-0.0002867	13	2572.729	-5119.310	98.66709	3.749143e-22

kable(Subset.dredged.extinct.model, "html", caption = "BAMM Extinction Male-bias Sexual Selection Dredge Table") %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

BAMM Extinction Male-bias Sexual Selection Dredge Table
	(Intercept)	bioclim4	log(range.size.m2)	PC1.LIG	residuals.PC1	Sexual_selection_ppca	bioclim4:Sexual_selection_ppca	log(range.size.m2):Sexual_selection_ppca	NPP:Sexual_selection_ppca	PC1.LIG:Sexual_selection_ppca	residuals.PC1:Sexual_selection_ppca	df	logLik	AICc	delta	weight
0	-2.835512	-7.3e-06	0.0004994	-0.0003236	0.0004112	-0.0003583	NA	NA	NA	NA	NA	8	2617.018	-5217.977	0.00000	9.982022e-01
8	-2.837029	-7.4e-06	0.0005575	-0.0003517	0.0004687	-0.0004496	NA	NA	NA	0.0004539	NA	9	2610.740	-5203.407	14.56967	6.846268e-04
16	-2.835892	-7.3e-06	0.0005110	-0.0003416	0.0003737	-0.0004570	NA	NA	NA	NA	-0.0002267	9	2610.652	-5203.231	14.74603	6.268432e-04
2	-2.835374	-7.2e-06	0.0004942	-0.0003057	0.0003976	-0.0048344	NA	0.0001612	NA	NA	NA	9	2610.384	-5202.694	15.28266	4.793257e-04
1	-2.835131	-7.0e-06	0.0004827	-0.0003044	0.0004030	-0.0008253	1.5e-06	NA	NA	NA	NA	9	2605.952	-5193.830	24.14720	5.697974e-06
24	-2.837371	-7.4e-06	0.0005678	-0.0003683	0.0004333	-0.0005410	NA	NA	NA	0.0004501	-0.0002114	10	2604.367	-5188.644	29.33283	4.262599e-07
10	-2.836889	-7.3e-06	0.0005522	-0.0003347	0.0004556	-0.0046804	NA	0.0001523	NA	0.0004513	NA	10	2604.101	-5188.112	29.86480	3.267087e-07
18	-2.835724	-7.3e-06	0.0005051	-0.0003243	0.0003679	-0.0041362	NA	0.0001330	NA	NA	-0.0001943	10	2604.018	-5187.947	30.03002	3.008033e-07
4	-2.835551	-7.3e-06	0.0005019	-0.0003257	0.0004112	-0.0001576	NA	NA	0	NA	NA	9	2602.884	-5187.695	30.28165	2.652418e-07
17	-2.835405	-6.9e-06	0.0004876	-0.0003165	0.0003449	-0.0012820	2.5e-06	NA	NA	NA	-0.0003181	10	2599.694	-5179.298	38.67931	3.982116e-09
9	-2.836626	-7.1e-06	0.0005397	-0.0003311	0.0004602	-0.0009541	1.6e-06	NA	NA	0.0004565	NA	10	2599.679	-5179.269	38.70807	3.925262e-09
3	-2.835149	-7.1e-06	0.0004842	-0.0002961	0.0003945	-0.0044284	1.0e-06	0.0001357	NA	NA	NA	10	2599.353	-5178.616	39.36092	2.832071e-09
26	-2.837205	-7.4e-06	0.0005620	-0.0003517	0.0004276	-0.0040318	NA	0.0001262	NA	0.0004485	-0.0001808	11	2597.730	-5173.352	44.62511	2.037042e-10
12	-2.837067	-7.4e-06	0.0005600	-0.0003538	0.0004688	-0.0002498	NA	NA	0	0.0004539	NA	10	2596.607	-5173.124	44.85298	1.817677e-10
20	-2.835901	-7.3e-06	0.0005117	-0.0003421	0.0003741	-0.0003955	NA	NA	0	NA	-0.0002247	10	2596.527	-5172.964	45.01281	1.678074e-10
6	-2.835390	-7.3e-06	0.0004952	-0.0003067	0.0003978	-0.0047065	NA	0.0001593	0	NA	NA	10	2596.258	-5172.427	45.55036	1.282580e-10
25	-2.836876	-6.9e-06	0.0005439	-0.0003425	0.0004039	-0.0013918	2.6e-06	NA	NA	0.0004525	-0.0003056	11	2593.413	-5164.719	53.25808	2.718774e-12
19	-2.835390	-6.9e-06	0.0004878	-0.0003117	0.0003459	-0.0028671	2.2e-06	0.0000611	NA	NA	-0.0002915	11	2593.127	-5164.146	53.83075	2.041833e-12
11	-2.836632	-7.1e-06	0.0005408	-0.0003235	0.0004522	-0.0042024	1.1e-06	0.0001223	NA	0.0004536	NA	11	2593.075	-5164.042	53.93520	1.937935e-12
5	-2.835167	-7.0e-06	0.0004849	-0.0003063	0.0004032	-0.0006555	1.5e-06	NA	0	NA	NA	10	2591.819	-5163.548	54.42891	1.514015e-12
28	-2.837380	-7.4e-06	0.0005686	-0.0003688	0.0004337	-0.0004705	NA	NA	0	0.0004502	-0.0002091	11	2590.242	-5158.376	59.60107	1.140276e-13
14	-2.836906	-7.3e-06	0.0005533	-0.0003357	0.0004558	-0.0045421	NA	0.0001503	0	0.0004513	NA	11	2589.975	-5157.843	60.13398	8.735578e-14
22	-2.835720	-7.3e-06	0.0005048	-0.0003240	0.0003677	-0.0041769	NA	0.0001335	0	NA	-0.0001950	11	2589.900	-5157.692	60.28495	8.100421e-14
27	-2.836861	-7.0e-06	0.0005441	-0.0003386	0.0004046	-0.0026825	2.3e-06	0.0000498	NA	0.0004516	-0.0002839	12	2586.845	-5149.562	68.41474	1.390423e-15
21	-2.835393	-6.9e-06	0.0004868	-0.0003158	0.0003443	-0.0013506	2.5e-06	NA	0	NA	-0.0003208	11	2585.575	-5149.042	68.93502	1.071937e-15
13	-2.836661	-7.1e-06	0.0005419	-0.0003330	0.0004604	-0.0007883	1.6e-06	NA	0	0.0004564	NA	11	2585.547	-5148.986	68.99156	1.042059e-15
7	-2.835165	-7.1e-06	0.0004852	-0.0002970	0.0003947	-0.0043083	1.0e-06	0.0001339	0	NA	NA	11	2585.227	-5148.347	69.62984	7.573414e-16
30	-2.837203	-7.4e-06	0.0005618	-0.0003516	0.0004275	-0.0040513	NA	0.0001264	0	0.0004485	-0.0001811	12	2583.611	-5143.095	74.88169	5.481093e-17
29	-2.836865	-6.9e-06	0.0005432	-0.0003419	0.0004034	-0.0014547	2.6e-06	NA	0	0.0004525	-0.0003080	12	2579.294	-5134.462	83.51548	7.312432e-19
23	-2.835374	-6.9e-06	0.0004867	-0.0003108	0.0003452	-0.0029967	2.2e-06	0.0000624	0	NA	-0.0002946	12	2579.012	-5133.896	84.08130	5.510563e-19
15	-2.836649	-7.1e-06	0.0005418	-0.0003245	0.0004524	-0.0040732	1.1e-06	0.0001205	0	0.0004536	NA	12	2578.949	-5133.771	84.20563	5.178431e-19
31	-2.836846	-7.0e-06	0.0005430	-0.0003377	0.0004039	-0.0027978	2.3e-06	0.0000509	0	0.0004516	-0.0002867	13	2572.729	-5119.310	98.66709	3.749143e-22

#Run model for DR
Subset.MCC.DR.top <- gls(MCC.DR ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = FALSE), 
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.DR.top, 'data/Subset.MCC.DR.top.rds')

#Run model for ND
Subset.MCC.ND.top <- gls(MCC.ND ~ Sexual_selection_ppca
                         + log(range.size.m2)
                         + bioclim4 #Seasonal variation
                         + residuals.PC1 #Spatial variation
                         + PC1.LIG #Long-term climate variation
                         + NPP,
                correlation = corPagel(1, phy = pruned.MCC.Subset.tree, fixed = TRUE), #lambda = 1
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.ND.top, 'data/Subset.MCC.ND.top.rds')

#Run the 100 models for DR and ND using the best model:
Subset.data.noMCC <- SS.subset %>% dplyr::select(TipLabel, Sexual_selection_ppca, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP)

#Take the restricted data and make it simpler with just responses and predictors.Note that we join the es.values for the 100 trees
Subset.DR.model.data <- lapply(es.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Subset.data.noMCC, 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "DR")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Subset.DR.model.data <- lapply(Subset.DR.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Prune the trees
Subset.pruned.trees<-lapply(passerine.trees, function(x) {
  drop.tip(x,x$tip.label[-match(SS.subset$TipLabel, x$tip.label)])
})

#Use mapply to create a list of PGLS global models
Subset.DR.pgls.models <- mcmapply(function(x,y) {
  gls(DR ~ Sexual_selection_ppca 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    correlation = corPagel(0.9711, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = Subset.DR.model.data, y = Subset.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.DR.pgls.models, "data/Subset.DR.pgls.models.rds")

#Now for Node Density:
Subset.ND.model.data <- lapply(nd.list, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Subset.data.noMCC, 
            x %>% as.data.frame() %>% tibble::rownames_to_column() %>% `colnames<-`(c("TipLabel", "ND")), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Subset.ND.model.data <- lapply(Subset.ND.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Subset.ND.pgls.models <- mcmapply(function(x,y) {
gls(ND ~ Sexual_selection_ppca 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    corPagel(1, phy = y, fixed = TRUE), 
    data = x, 
    method = "REML")
}, x = Subset.ND.model.data, y = Subset.pruned.trees, 
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.ND.pgls.models, "data/Subset.ND.pgls.models.rds")

#BAMM Top Models 

Subset.MCC.Lambda.top <- gls(log(mean.lambda) ~ Sexual_selection_ppca
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ I(1/log(var.lambda)),
                correlation = corBrownian(phy = pruned.MCC.Subset.tree), 
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.Lambda.top, 'data/Subset.MCC.Lambda.top.rds')

Subset.MCC.Mu.top <- gls(log(mean.mu) ~ Sexual_selection_ppca
                      + log(range.size.m2)
                      + bioclim4 #Seasonal variation
                      + residuals.PC1 #Spatial variation
                      + PC1.LIG #Long-term climate variation
                      + NPP,
                weights = ~ I(1/log(var.mu)),
                correlation = corBrownian(phy = pruned.MCC.Subset.tree), 
                data = SS.subset, 
                method = "REML")
saveRDS(Subset.MCC.Mu.top, 'data/Subset.MCC.Mu.top.rds')

#Now for BAMM 100 models

Subset.BAMM.model.data <- lapply(BAMM.df, function(x) { #es.list is a list of ES values calculated earlier
  left_join(Subset.data.noMCC %>% dplyr::select(TipLabel, Sexual_selection_ppca, range.size.m2, bioclim4, residuals.PC1, PC1.LIG, NPP), 
            x %>% as.data.frame(), 
            by = "TipLabel")
})

#PGLS needs tiplabel as rowname
Subset.BAMM.model.data <- lapply(Subset.BAMM.model.data, function(x) {
  tibble::column_to_rownames(x, "TipLabel")})

#Use mapply to create a list of PGLS global models
Subset.BAMM.lambda.pgls.models <- mcmapply(function(x,y) {
  gls(log(mean.lambda) ~ Sexual_selection_ppca 
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ I(1/log(var.lambda)),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Subset.BAMM.model.data, y = Subset.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.BAMM.lambda.pgls.models, "data/Subset.BAMM.lambda.pgls.models.rds")

#Use mapply to create a list of PGLS global models
Subset.BAMM.mu.pgls.models <- mcmapply(function(x,y) {
  gls(log(mean.mu) ~ Sexual_selection_ppca
         + log(range.size.m2)
         + bioclim4 #Seasonal variation
         + residuals.PC1 #Spatial variation
         + PC1.LIG #Long-term climate variation
         + NPP,
    weights = ~ I(1/log(var.mu)),      
    corBrownian(phy = y), #lambda = 1
    data = x, 
    method = "REML")
}, x = Subset.BAMM.model.data, y = Subset.pruned.trees,
SIMPLIFY = FALSE, #Prevents the catastrophic loss of structure, names and attributes
mc.cores = 8) #Specify core number 

saveRDS(Subset.BAMM.mu.pgls.models, "data/Subset.BAMM.mu.pgls.models.rds")

The phylogenetic signal from the DR model is ~0.97. The other models had phylogenetic signals ~ 1, as such they were fixed at 1 to avoid problems of convergence.

Subset.MCC.DR.top <- readRDS('data/Subset.MCC.DR.top.rds')
Subset.MCC.DR.top[["modelStruct"]][["corStruct"]] %>% `names<-`("DR lambda") %>% pander()

DR lambda: 0.9711

Subset.DR.pgls.models <- readRDS("data/Subset.DR.pgls.models.rds")
Subset.MCC.DR.top <- readRDS('data/Subset.MCC.DR.top.rds')

Subset.cols <- brewer.pal(n = 7, name = "Dark2")[-6]

Subset.DR.pgls.summary <- lapply(Subset.DR.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.MCC.DR.summary <- data.frame(Subset.MCC.DR.top$coefficients, confint(Subset.MCC.DR.top)) %>% tibble::rownames_to_column()

Subset.DR.pgls.summary <- bind_rows(Subset.DR.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Subset.MCC.DR.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

parameter_names <- c(
                    `bioclim4` = "Temperature Seasonality",
                    `log(range.size.m2)` = "Range Size (log-transformed)",
                    `NPP` = "NPP",
                    `PC1.LIG` = "Long-term Temperature Variation",
                    `residuals.PC1` = "Spatial Temperature Variation",
                    `Sexual_selection_ppca` = "Sexual Selection"
                    )

Subset.DR.pgls.summary$Parameter = factor(Subset.DR.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.DR.summary$Parameter = factor(Subset.MCC.DR.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.DR.plot <-Subset.DR.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.DR.plot <- Subset.DR.plot + geom_errorbarh(data = Subset.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.DR.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS DR")

####################################
#For ND

Subset.ND.pgls.models <- readRDS("data/Subset.ND.pgls.models.rds")
Subset.MCC.ND.top <- readRDS('data/Subset.MCC.ND.top.rds')

Subset.ND.pgls.summary <- lapply(Subset.ND.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.MCC.ND.summary <- data.frame(Subset.MCC.ND.top$coefficients, confint(Subset.MCC.ND.top)) %>% tibble::rownames_to_column()

Subset.ND.pgls.summary <- bind_rows(Subset.ND.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Subset.MCC.ND.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Subset.ND.pgls.summary$Parameter = factor(Subset.ND.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.ND.summary$Parameter = factor(Subset.MCC.ND.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.ND.plot <-Subset.ND.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.ND.plot <- Subset.ND.plot + geom_errorbarh(data = Subset.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.ND.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS ND")

######################

#For Lambda
Subset.BAMM.lambda.pgls.models <- readRDS("data/Subset.BAMM.lambda.pgls.models.rds")
Subset.MCC.Lambda.top <- readRDS('data/Subset.MCC.Lambda.top.rds')

Subset.Lambda.pgls.summary <- lapply(Subset.BAMM.lambda.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.MCC.Lambda.summary <- data.frame(Subset.MCC.Lambda.top$coefficients, confint(Subset.MCC.Lambda.top)) %>% tibble::rownames_to_column()

Subset.Lambda.pgls.summary <- bind_rows(Subset.Lambda.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Subset.MCC.Lambda.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Subset.Lambda.pgls.summary$Parameter = factor(Subset.Lambda.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.Lambda.summary$Parameter = factor(Subset.MCC.Lambda.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.Lambda.plot <-Subset.Lambda.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.Lambda.plot <- Subset.Lambda.plot + geom_errorbarh(data = Subset.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.Lambda.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS Speciaton")

######################

#For Mu

Subset.BAMM.Mu.pgls.models <- readRDS("data/Subset.BAMM.mu.pgls.models.rds")
Subset.MCC.Mu.top <- readRDS('data/Subset.MCC.Mu.top.rds')

Subset.Mu.pgls.summary <- lapply(Subset.BAMM.Mu.pgls.models, function(x) {
  data.frame(x$coefficients, confint(x)) %>% tibble::rownames_to_column()
})

Subset.MCC.Mu.summary <- data.frame(Subset.MCC.Mu.top$coefficients, confint(Subset.MCC.Mu.top)) %>% tibble::rownames_to_column()

Subset.Mu.pgls.summary <- bind_rows(Subset.Mu.pgls.summary) %>% `colnames<-`(c("Parameter", "Estimate", "LCI", "UCI"))

colnames(Subset.MCC.Mu.summary) <- c("Parameter", "Estimate", "LCI", "UCI")

Subset.Mu.pgls.summary$Parameter = factor(Subset.Mu.pgls.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.MCC.Mu.summary$Parameter = factor(Subset.MCC.Mu.summary$Parameter, levels=c('bioclim4','log(range.size.m2)','NPP','PC1.LIG', 'residuals.PC1', 'Sexual_selection_ppca'))

Subset.Mu.plot <-Subset.Mu.pgls.summary %>% filter(Parameter != "(Intercept)") %>% 
  ggplot(aes(x = Estimate, y = Parameter, fill = Parameter)) +
  geom_density(aes(y = ..scaled..),
               position = position_nudge(y= 0.6))+
  geom_jitter(shape = 21, alpha = 0.5, size = 0.75)+
  geom_vline(xintercept = 0)+
  theme_minimal()+
  theme(axis.text.y = element_blank(),
        legend.position = "none")

Subset.Mu.plot <- Subset.Mu.plot + geom_errorbarh(data = Subset.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
                 aes(xmin = LCI,
                     xmax = UCI, y = Parameter,
                     color = Parameter), 
                 height = 0, show.legend = F, 
                 position = position_nudge(y= -0.75))+
  
  geom_point(data = Subset.MCC.Mu.summary %>% filter(Parameter != "(Intercept)"), 
             aes(x = Estimate, y = Parameter, fill = Parameter), 
             shape=21, color = "grey20",
             size = 3,
             position = position_nudge(y= -0.75)) +
  scale_y_discrete(expand = c(0.5,.5))+
  facet_wrap(~ Parameter, scales = "free", nrow = 6, labeller = as_labeller(parameter_names))+
  scale_fill_manual(values = Subset.cols)+
  scale_color_manual(values = Subset.cols)+
  theme(panel.grid.major.y = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 7),
        plot.margin=unit(c(.5,.5,.5,.5),"cm"))+
  ggtitle("Male-Bias SS Extinction")

#########################################

grid.arrange(symmetrise_scale(Subset.DR.plot, "x"),
             symmetrise_scale(Subset.ND.plot, "x"),
             symmetrise_scale(Subset.Lambda.plot, "x"), 
             symmetrise_scale(Subset.Mu.plot, "x"), 
             nrow = 1)

Figure S11: Measures of sexual selection a phylogenetic pca of sexual dimorphism, social polygyny and paternal care provide some evidence to suggest sexual selection predicts extinction rate (DR and ND models). This figure is seen in Figure 1 within the manuscript

Table S11: Highest posterior density intervals are calculated for the above figure (ie; models using Male-bias sexual selection measures). using this measure the models using \(\lambda_{DR}\) have model estimates that fall on the positive side 95 % of the time in 100 trees. For \(\lambda_{ND}\), there is a positive skew, however the 95 % HPD interval overlaps zero. \(\lambda_{BAMM}\) also shows a lesser positive skew. The HPD range is determined using the hdi function with a 95 % credible interval. These intervals do not take into account the variance associated with each interval and thus are not an estimate of model precision. Intervals not overlapping zero suggest that 95 % of trees from the posterior generate a model estimate for the given parameter that in the same direction (+ or -). These intervals are calculated in the same way as in Table S7 and Table S9.

Subset.hpd.DR.top <- list()
Subset.DR.pgls.summary <- na.omit(Subset.DR.pgls.summary)
for(x in unique(Subset.DR.pgls.summary$Parameter)) {
Subset.hpd.DR.top[[x]] = hdi(Subset.DR.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.DR.top <- bind_rows(Subset.hpd.DR.top) %>% `rownames<-`(c("Lower", "Upper"))


saveRDS(Subset.hpd.DR.top, 'data/Subset.hpd.DR.top.rds')

#For ND
Subset.hpd.ND.top <- list()
Subset.ND.pgls.summary <- na.omit(Subset.ND.pgls.summary)
for(x in unique(Subset.ND.pgls.summary$Parameter)) {
Subset.hpd.ND.top[[x]] = hdi(Subset.ND.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.ND.top <- bind_rows(Subset.hpd.ND.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Subset.hpd.ND.top, 'data/Subset.hpd.ND.top.rds')

###############
Subset.hpd.Lambda.top <- list()
Subset.Lambda.pgls.summary <- na.omit(Subset.Lambda.pgls.summary)
for(x in unique(Subset.Lambda.pgls.summary$Parameter)) {
Subset.hpd.Lambda.top[[x]] = hdi(Subset.Lambda.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.Lambda.top <- bind_rows(Subset.hpd.Lambda.top) %>% `rownames<-`(c("Lower", "Upper")) 
saveRDS(Subset.hpd.Lambda.top, 'data/Subset.hpd.Lambda.top.rds')

################
Subset.hpd.Mu.top <- list()
Subset.Mu.pgls.summary <- na.omit(Subset.Mu.pgls.summary)
for(x in unique(Subset.Mu.pgls.summary$Parameter)) {
Subset.hpd.Mu.top[[x]] = hdi(Subset.Mu.pgls.summary %>% filter(Parameter == x) %>% dplyr::select(Estimate))
}
Subset.hpd.Mu.top <- bind_rows(Subset.hpd.Mu.top) %>% `rownames<-`(c("Lower", "Upper"))

saveRDS(Subset.hpd.Mu.top, 'data/Subset.hpd.Mu.top.rds')

Subset.hpd.DR.top %>% pander(split.table = Inf, digits = 3, caption = "Subset DR HPD Interval")

Subset DR HPD Interval
	Sexual_selection_ppca	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	0.00911	-0.0158	-9.31e-05	0.000313	-0.0036	-2.96e-06
Upper	0.0608	-0.000501	2.76e-05	0.0128	0.00894	1.19e-06

Subset.hpd.ND.top %>% pander(split.table = Inf, digits = 3, caption = "Subset ND HPD Interval")

Subset ND HPD Interval
	Sexual_selection_ppca	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.00031	-0.000426	-3.63e-06	-3.05e-05	-0.000249	-9.23e-08
Upper	0.00154	0.000166	1.64e-06	0.000487	0.000249	6.61e-08

Subset.hpd.Lambda.top %>% pander(split.table = Inf, digits = 3, caption = "Subset BAMM Speciation HPD Interval")

Subset BAMM Speciation HPD Interval
	Sexual_selection_ppca	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.0198	-0.00853	-7.27e-05	-0.0114	-0.0106	-1.23e-06
Upper	0.0313	0.00801	5.36e-05	0.00827	0.00395	1.91e-06

Subset.hpd.Mu.top %>% pander(split.table = Inf, digits = 3, caption = "Subset BAMM Extinction HPD Interval")

Subset BAMM Extinction HPD Interval
	Sexual_selection_ppca	log(range.size.m2)	bioclim4	residuals.PC1	PC1.LIG	NPP
Lower	-0.0441	-0.0134	-0.000145	-0.0125	-0.0156	-1.97e-06
Upper	0.0431	0.013	7.66e-05	0.0104	0.00951	3.36e-06

Additional Figures

MCC.BAMM.model.data <- readRDS('data/MCC.BAMM.model.data.rds')
restricted.data %>% ggplot(aes(x = MCC.DR, y = MCC.BAMM.model.data$mean.lambda))+
  geom_point()+
  scale_y_continuous(trans = "log")+
  theme_minimal()

R Session information

sessionInfo() %>% pander

R version 3.3.1 (2016-06-21)

**Platform:** x86_64-apple-darwin13.4.0 (64-bit)

locale: en_AU.UTF-8||en_AU.UTF-8||en_AU.UTF-8||C||en_AU.UTF-8||en_AU.UTF-8

attached base packages: parallel, grid, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: bindrcpp(v.0.2.2), ggtree(v.1.6.11), EBImage(v.4.16.0), RColorBrewer(v.1.1-2), kableExtra(v.0.7.0), phangorn(v.2.4.0), HDInterval(v.0.2.0), MuMIn(v.1.42.1), coda(v.0.19-1), car(v.3.0-2), carData(v.3.0-1), BAMMtools(v.2.1.6), dplyr(v.0.7.6), ggExtra(v.0.8), reshape2(v.1.4.3), purrr(v.0.2.4), caper(v.1.0.1), mvtnorm(v.1.0-6), MASS(v.7.3-50), ggridges(v.0.5.0), brms(v.2.4.0), Rcpp(v.0.12.18), phytools(v.0.6-44), maps(v.3.2.0), gridExtra(v.2.3), geiger(v.2.0.6), lme4(v.1.1-18-1), Matrix(v.1.2-14), repmis(v.0.5), diversitree(v.0.9-10), ape(v.5.1), mgcv(v.1.8-24), nlme(v.3.1-128), tidyr(v.0.7.2), knitr(v.1.20), pander(v.0.6.2) and ggplot2(v.3.0.0)

loaded via a namespace (and not attached): readxl(v.1.0.0), backports(v.1.1.2), fastmatch(v.1.1-0), plyr(v.1.8.4), igraph(v.1.1.2), lazyeval(v.0.2.1), splines(v.3.3.1), crosstalk(v.1.0.0), rstantools(v.1.5.1), inline(v.0.3.15), digest(v.0.6.16), htmltools(v.0.3.6), tiff(v.0.1-5), rsconnect(v.0.8.8), gdata(v.2.18.0), magrittr(v.1.5), openxlsx(v.4.1.0), readr(v.1.1.1), matrixStats(v.0.54.0), R.utils(v.2.7.0), xts(v.0.11-0), jpeg(v.0.1-8), colorspace(v.1.3-2), rvest(v.0.3.2), haven(v.1.1.2), jsonlite(v.1.5), bindr(v.0.1.1), survival(v.2.39-4), zoo(v.1.8-3), glue(v.1.3.0), gtable(v.0.2.0), R.cache(v.0.13.0), rstan(v.2.17.3), BiocGenerics(v.0.20.0), abind(v.1.4-5), scales(v.1.0.0), miniUI(v.0.1.1.1), plotrix(v.3.7-3), viridisLite(v.0.3.0), xtable(v.1.8-2), foreign(v.0.8-66), subplex(v.1.4-1), deSolve(v.1.20), stats4(v.3.3.1), StanHeaders(v.2.17.2), DT(v.0.4), animation(v.2.5), htmlwidgets(v.1.2), httr(v.1.3.1), threejs(v.0.3.1), gplots(v.3.0.1), pkgconfig(v.2.0.2), loo(v.2.0.0), R.methodsS3(v.1.7.1), locfit(v.1.5-9.1), labeling(v.0.3), tidyselect(v.0.2.4), rlang(v.0.2.2), later(v.0.7.3), munsell(v.0.5.0), cellranger(v.1.1.0), tools(v.3.3.1), fftwtools(v.0.9-8), evaluate(v.0.10.1), stringr(v.1.2.0), yaml(v.2.2.0), zip(v.1.0.0), caTools(v.1.17.1), mime(v.0.5), R.oo(v.1.22.0), xml2(v.1.2.0), bayesplot(v.1.6.0), shinythemes(v.1.1.1), png(v.0.1-7), curl(v.3.2), clusterGeneration(v.1.3.4), tibble(v.1.3.4), stringi(v.1.1.6), Brobdingnag(v.1.2-6), forcats(v.0.3.0), lattice(v.0.20-33), nloptr(v.1.0.4), markdown(v.0.8), shinyjs(v.1.0), msm(v.1.6.6), combinat(v.0.0-8), bridgesampling(v.0.5-2), data.table(v.1.11.4), bitops(v.1.0-6), httpuv(v.1.4.5), R6(v.2.2.2), promises(v.1.0.1), KernSmooth(v.2.23-15), rio(v.0.5.10), codetools(v.0.2-14), colourpicker(v.1.0), gtools(v.3.8.1), assertthat(v.0.2.0), rprojroot(v.1.3-2), withr(v.2.1.2), shinystan(v.2.4.0), mnormt(v.1.5-5), expm(v.0.999-2), hms(v.0.4.2), quadprog(v.1.5-5), minqa(v.1.2.4), rmarkdown(v.1.10), numDeriv(v.2016.8-1), scatterplot3d(v.0.3-41), shiny(v.1.1.0), base64enc(v.0.1-3) and dygraphs(v.1.1.1.6)

References

Armenta, J. K., P. O. Dunn, and L. A. Whittingham. 2008. Quantifying avian sexual dichromatism: A comparison of methods. Journal of Experimental Biology 211:2423. Journal Article.

BirdLife International and Handbook of the Birds of the World. 2017. Bird species distribution maps of the world. http://datazone.birdlife.org/species/requestdis.

Dale, J., C. J. Dey, K. Delhey, B. Kempenaers, and M. Valcu. 2015. The effects of life history and sexual selection on male and female plumage colouration. Nature 527:367–370. Journal Article.

Del Hoyo, J., A. Elliott, and D. Christie. n.d. Handbook of the birds of the world (Vols. 8-16). Book, Lynx Edicions 2003-2011.

Harvey, M. G., G. F. Seeholzer, B. T. Smith, D. L. Rabosky, A. M. Cuervo, and R. T. Brumfield. 2017. Positive association between population genetic differentiation and speciation rates in new world birds. Proceedings of the National Academy of Sciences 114:6328–6333.

Harvey Michael, G., L. Rabosky Daniel, and N. Cooper. 2017. Continuous traits and speciation rates: Alternatives to state‐dependent diversification models. Methods in Ecology and Evolution 9:984–993. Journal Article.

Jetz, W., G. H. Thomas, J. B. Joy, K. Hartmann, and A. O. Mooers. 2012. The global diversity of birds in space and time. Nature 491:444–448. Journal Article.

Quintero, I., and W. Jetz. 2018. Global elevational diversity and diversification of birds. Nature 555:246.

Schliep, K. P. 2010. Phangorn: Phylogenetic analysis in r. Bioinformatics 27:592–593.

Sexual selection and speciation: a macro-evolutionary approach

Supplementary Material

Justin G. Cally¹, Iliana Medina Guzmán¹ Luke Holman¹, Devi Stuart-Fox¹

¹The University of Melbourne

Compiling the data

Speciation Rate

Sexual Selection

Environmental predictors

Generating biologically relevent predictors of extinction

Spatial Variability PCAs

Long term climate variability (LIG anomaly)

Correlations between environmental predictors

Sexual Dimorphism

Analysis

Simple speciation measures: Diversification rate (DR) and Node Density (ND)

BAMM measures of speciation and extinction

Set up BAMM parameters

Run analysis

Read in the event data and extract tip data

Analysis of BAMM results

Credible number of shifts

PGLS Models

PGLS Models on DR and ND

PGLS Models on BAMM estimates

PGLS Model Summary

Subsetted analysis with Armenta data

Analysis using male-bias measure of sexual selection

Additional Figures

R Session information

References

Sexual selection and speciation: a macro-evolutionary approach

Supplementary Material

Justin G. Cally1, Iliana Medina Guzmán1 Luke Holman1, Devi Stuart-Fox1 1The University of Melbourne

Compiling the data

Speciation Rate

Sexual Selection

Environmental predictors

Generating biologically relevent predictors of extinction

Spatial Variability PCAs

Long term climate variability (LIG anomaly)

Correlations between environmental predictors

Sexual Dimorphism

Analysis

Simple speciation measures: Diversification rate (DR) and Node Density (ND)

BAMM measures of speciation and extinction

Set up BAMM parameters

Run analysis

Read in the event data and extract tip data

Analysis of BAMM results

Credible number of shifts

PGLS Models

PGLS Models on DR and ND

PGLS Models on BAMM estimates

PGLS Model Summary

Subsetted analysis with Armenta data

Analysis using male-bias measure of sexual selection

Additional Figures

R Session information

References

Justin G. Cally¹, Iliana Medina Guzmán¹ Luke Holman¹, Devi Stuart-Fox¹

¹The University of Melbourne